Validating genomic offset predictions with mortality and height data from common gardens

Published

January 16, 2025

1 Introduction

In this report, we evaluate whether genomic offset (GO) predictions from the different GEA methods (LFMM, GF, GDM and RDA) and SNP sets (control, candidate and all SNPs for LFMM) are good predictors of fitness proxies (height and mortality data) from five clonal common gardens located in Spain (Asturias, Madrid, Cáceres), Portugal (Fundão) and France (Pierroton). For that, we predicted the genomic offset of each population when transplanted in the environment of the common garden (so instead of predicting the genomic offset into future climates, we predict it in the new environment of the common garden, i.e. space-for-time approach). The rational is that populations with higher genomic offset in a given common garden are expected to have lower relative fitness in the common garden.

We also calculate the climatic transfer distances (CTD) between the climate-of-origin of the populations and the climate of the common garden (i.e. the absolute value of the difference between the climate of origin of the populations and the climate in the common garden between the tree planting date and the measurement date). We then estimate the association between tree height and mortality and the climatic transfer distances. The climatic transfer distances are calculated for the same set of climatic variables as the one used to calculate the genomic offset (one CTD per climatic variable).

2 Data

2.1 Genomic offset predictions

We load the genomic offset predictions estimated from the different GEAs.

Code
# Codes of the SNP sets
snp_sets <- readRDS(here("outputs/list_snp_sets.rds"))

# Which RDA method?
RDA_method <- "predict"

# List of methods
meth_names <- c("GDM", "GF", "LFMM", "pRDA", "RDA")

df <- lapply(snp_sets, function(x){

list_snps <- list()

list_snps$GDM <- readRDS(file=here("outputs/GDM/go_predictions.rds"))[[x$set_code]][["go_cg"]]
list_snps$GF <- readRDS(file=here("outputs/GF/go_predictions.rds"))[[x$set_code]][["go_cg"]]
list_snps$LFMM <- readRDS(file=here("outputs/LFMM/go_predictions.rds"))[[x$set_code]][["go_cg"]]
list_snps$pRDA <- readRDS(file=here(paste0("outputs/RDA/go_predictions_",RDA_method,".rds")))[[x$set_code]][["go_cg_corrected"]]
list_snps$RDA <- readRDS(file=here(paste0("outputs/RDA/go_predictions_",RDA_method,".rds")))[[x$set_code]][["go_cg_uncorrected"]]


list_snps %>% 
  bind_rows(.id="method") %>% 
  mutate(input_code = x$set_code,
         input_name = x$set_name,
         method_input_code = paste0(method, "_", input_code),
         method_input_name = paste0(method, " - ", input_name))

}) %>% 
  bind_rows() %>% 
  pivot_longer(cols=names(cg_names),names_to="cg",values_to="varX") 

#####


# df <- lapply(c("control","cand"), function(x){
# 
# list_snps <- list()
# 
# list_snps[[x]]$GDM <- readRDS(file=here("outputs/GDM/go_predictions.rds"))[[x]][["go_cg"]]
# list_snps[[x]]$GF <- readRDS(file=here("outputs/GF/go_predictions.rds"))[[x]][["go_cg"]]
# list_snps[[x]]$RDA <- readRDS(file=here("outputs/RDA/go_predictions.rds"))[[x]][["go_cg"]]
# list_snps[[x]]$LFMM <- readRDS(file=here("outputs/LFMM/go_predictions_snpsets.rds"))[[x]][["go_cg"]]
# 
# 
# list_snps <- list_snps[[x]] %>% 
#   bind_rows(.id="method_type") %>% 
#   mutate(method_input = x)
# 
# return(list_snps)
# }) %>% bind_rows()
# 
# 
# df <- readRDS(file=here("outputs/LFMM/go_predictions_allsnps.rds"))[["go_cg"]] %>% 
#    mutate(method_type = "LFMM",
#           method_input = "all") %>% 
#   bind_rows(df) %>% 
#   pivot_longer(cols=c(readRDS(here("data/ClimaticData/CommonGardens/ClimateCG.rds"))[["cg"]]),names_to="cg",values_to="varX") %>% 
#   mutate(method = paste0(method_type, "_",method_input))

2.2 Climatic transfer distances

We calculate the climatic transfer distances.

Code
# selected climatic variables
clim_var <- readRDS(here("data/ClimaticData/NamesSelectedVariables.rds"))

# climatic data in the common gardens (btw planting and measurement dates)
cg_clim <-readRDS(here("data/ClimaticData/CommonGardens/ClimateCG.rds")) %>% 
  dplyr::select(cg,any_of(clim_var)) 

# Loading point estimate climatic data at the location of the populations (reference climatic period)
adj <- "noADJ"  # not adjusted for elevation
ref_period <- "ref_1901_1950" # reference period 1901-1950
clim_past <- readRDS(here(paste0("data/ClimaticData/MaritimePinePops/ClimatePopulationLocationPointEstimates_ReferencePeriods_",adj,".rds")))[[ref_period]]$ref_means %>%
  dplyr::select(pop,any_of(clim_var))

# Preparing the dataset for 
clim_dist <- clim_past

df <- lapply(cg_clim[["cg"]], function(x){
  
# Calculating the CTD for each climatic variable
for(var in clim_var){clim_dist[[var]] <- ( clim_past[[var]] - cg_clim %>% filter(cg == x) %>%  pull(var) ) %>% abs()} 
#names(clim_dist) <- c("pop", str_c(names(clim_past[,-1]),"_CTD"))

# Calculating the Euclidean distances for all climatic variables
list_scaled_clim_df <- generate_scaled_clim_datasets(clim_var, clim_pred = cg_clim)
cg_clim_i <- list_scaled_clim_df$clim_pred %>% filter(cg == x) %>% dplyr::select(-cg)

clim_dist$EucliDist <- map_dbl(1:nrow(list_scaled_clim_df$clim_ref[,-1]),
                               ~ sqrt(sum((list_scaled_clim_df$clim_ref[.x,-1] - cg_clim_i)^2)))
   
return(clim_dist)
}) %>%  
  setNames(cg_clim[["cg"]]) %>% 
  bind_rows(.id="cg") %>% 
  pivot_longer(cols=c(any_of(clim_var),"EucliDist"),names_to="input_code",values_to="varX") %>% 
  mutate(method = "CTD",
         method_input_code = paste0("CTD_",input_code),
         input_name = case_when(input_code == "bio1" ~ "Mean annual temperature (bio1, °C)",
                                input_code == "bio12" ~ "Annual precipitation (bio12, mm)",
                                input_code == "bio15" ~ "Precipitation seasonality (bio15, index)",
                                input_code == "bio3" ~ "Isothermality (bio3, index)",
                                input_code == "bio4" ~ "Temperature seasonality (bio4, °C)",
                                input_code == "SHM" ~ "Summer heat moisture index (SHM, °C/mm)",
                                input_code == "EucliDist" ~ "All climatic variables (Euclidean distance)"),
         method_input_code = paste0("CTD_", input_code),
         method_input_name = paste0("CTD - ", input_name)) %>% 
  bind_rows(df)

2.3 Correlation CTDs vs GO predictions

We calculate the correlations among GO predictions and climatic transfer distances in the common gardens.

Code
lapply(unique(df$cg), function(site_i){
  
df_corr <- df %>% 
  dplyr::filter(cg==site_i) %>% 
  dplyr::select(method_input_code, varX,pop) %>% 
  pivot_wider(values_from = varX, names_from = method_input_code)

# Corrplot with ggplot
# ====================

correlations <- df_corr %>%  dplyr::select(-pop) %>% cor

# Order variables alphabetically
correlations <- correlations[order(rownames(correlations)), order(colnames(correlations))]
    
# Melt the correlation matrix into a long format
cor_data <- reshape2::melt(correlations)
  
# Filter to only show the lower triangle
cor_data <- cor_data %>% filter(as.numeric(Var1) > as.numeric(Var2))
  
# Plot
ggplot(cor_data, 
       aes(Var1, Var2, fill = value)) +
  geom_tile(color = "white") +
  scale_fill_gradient2(low = "#BB4444", mid = "#FFFFFF", high = "#4477AA", midpoint = 0, limits = c(-1, 1), name="Correlation") +
  geom_text(aes(label = round(value, 2)), color = "black", size = 3) +
  theme_minimal() +
  theme(axis.text.x = element_text(size=11,
                                   angle = 45, 
                                   hjust = 1),
        axis.text.y = element_text(size=11)) +
  labs(x = NULL, y = NULL)
  
# Corrplot with make_corrplot_pdf
# ===============================
fig_options <- list(
  path = here(paste0("figs/ValidationCommonGarden/Corrplot_GOpreds_CTDs_",site_i,".pdf")),
  width=10,
  height=8)

make_corrplot_pdf(df = df_corr,
                  variables = colnames(df_corr)[-c(1:2)],
                  text_size = 0.5,
                  number_size = 0.4,
                  fig_options = fig_options)
  
})

2.4 Phenotypic data

We load the survival and mortality data from the five common gardens.

Code
pheno_data <- readRDS(file=here("data/CommonGardenData/PhenoDataNovember2019_AnnualTraits_UpdatedSept2021_AllSites.rds")) %>% dplyr::rename(pop=prov)

no_nas <- sapply(pheno_data, function(x) length(x)-sum(is.na(x)))

list_pheno <- list()

list_pheno$`Asturias, Spain (37 months)` <- table(pheno_data$site,pheno_data$AST_survmar14)["asturias",]
list_pheno$`Bordeaux, France (85 months)` <- table(pheno_data$site,pheno_data$BDX_surv18)["bordeaux",]
list_pheno$`Cáceres, Spain (8 months)` <- table(pheno_data$site,pheno_data$CAC_survdec11)["caceres",]
list_pheno$`Madrid, Spain (13 months)` <- table(pheno_data$site,pheno_data$MAD_survdec11)["madrid",]
list_pheno$`Fundão, Portugal (27 months)` <- table(pheno_data$site,pheno_data$POR_survmay13)["portugal",]

df_exp_design <- list_pheno %>% 
  bind_rows(.id="cg") %>% 
  setNames(c("Common garden (tree age)",
             "Number of dead trees",
             "Number of trees alive")) %>% 
  mutate("Number of height measurements"=c(no_nas[["AST_htmar14"]],
                                           no_nas[["BDX_htnov18"]],
                                           no_nas[["CAC_htdec11"]],
                                           no_nas[["MAD_htdec11"]],
                                           no_nas[["POR_htmay13"]]))

saveRDS(df_exp_design, file = here("tables/ExpDesign_CG.rds"))

df_exp_design  %>% kable_mydf
Common garden (tree age) Number of dead trees Number of trees alive Number of height measurements
Asturias, Spain (37 months) 206 3976 3976
Bordeaux, France (85 months) 119 3225 3217
Cáceres, Spain (8 months) 3818 366 340
Madrid, Spain (13 months) 3138 1046 1046
Fundão, Portugal (27 months) 1517 2666 2665

Height and survival measurements in each common garden:

  • Asturias (Spain) in March 2014 when the trees were 37 month-old (trees were planted in February 2011).

  • Bordeaux (France) in November 2018 when the trees were 85 month-old (trees were planted in October 2011).

  • Cáceres (Spain) in December 2011 when the trees were 8 month-old (trees were planted in April 2011). Note that for this common garden, we calculate the bioclimatic variables for the entire year 2011 (instead of calculating the variables only for months between April and December). Indeed, the calculation of the annual bioclimatic variables will be wrong if we do not account for some months, e.g. the mean annual temperature will be higher than expected because we do not account for some winter months.

  • Madrid (Spain) in December 2011 when the trees were 13 month-old (trees were planted in November 2010).

  • Fundão (Portugal) in May 2013 when the trees were 27 month-old (trees were planted in February 2011)

3 Mortality models

In this section, we want to determine whether genomic offset (GO) or climate transfer distances (CTD) are associated with the proportion of dead trees in the populations, independently in the five common gardens. We build a model that assumes that the initial sapling height acts as a confounder. Indeed, trees that were higher at the time of planting have a higher probability of survival. This is particularly true in Madrid and Cáceres (Spain) where the seedlings experienced an extreme drought event the same year they were planted, resulting in very high mortality rates. Here is the model:

\[\begin{align*} a_{p} &\sim \text{Binomial} (N_{p},p_{p}) \\ \text{logit}(p_{p}) &= \beta_{0} + \beta_{H}H_{p} + \beta_{X}X_{p}\\ \end{align*}\]

with \(a_{p}\) the count of individuals that died in the population \(p\), \(N_{p}\) the total number of individuals in the population \(p\) (=number of individuals that were initially planted in the common garden), \(p_p\) is the estimated probability of mortality in the population \(p\), \(X_{p}\) is the genomic offset or climatic transfer distance for the population \(p\) and \(H_{p}\) is the BLUPs for height of the population \(p\) (population varying intercepts calculated across all common gardens in the model 1 of Archambeau et al. 2022). \(H_{p}\) is used as a proxy of the initial tree height.

We used the following weakly informative priors:

\[\begin{align*} \begin{bmatrix} \beta_{0,c} \\ \beta_{H} \\ \beta_{X} \end{bmatrix} &\sim \mathcal{N}(0,5) \end{align*}\]

3.1 The variables

3.1.1 Initial tree height (confounder)

We load the population-specific intercepts from the model 1 of Archambeau et al. (2022).

Code
mod_arch2022 <- readRDS(file=here("data/Archambeauetal2022_MOD1.rds"))
pop_heights <- mod_arch2022$fit %>% 
  broom.mixed::tidyMCMC(estimate.method = "mean",conf.int = T) %>% # we take the mean of the prov random intercepts
  filter(str_detect(term, "^(r_prov\\[)")) %>% 
  dplyr::rename(height=estimate,pop=term) %>% 
  mutate(pop=str_sub(pop,8,-12))

pop_heights[1:5,] %>% kable_mydf
pop height std.error conf.low conf.high
ALT 0.09 0.04 0.01 0.17
ARM 0.12 0.04 0.04 0.20
ARN -0.09 0.03 -0.16 -0.03
BAY -0.14 0.03 -0.20 -0.07
BON -0.04 0.04 -0.12 0.04

DRYAD REPOSITORY: We export the height intercepts from Archambeau et al. (2022) to the DRYAD repository: HeightIntercepts_Archambeauetal2022.csv.

Code
pop_heights %>% write_csv(here("data/DryadRepo/HeightIntercepts_Archambeauetal2022.csv"))

3.1.2 Survival data

The response variable is counts of dead trees. To calculate these counts, we load the survival data in the common gardens, in which 0 corresponds to dead trees and 1 to survivors.

Code
surv_measurements <- c("AST_survmar14","BDX_surv18","CAC_survdec11","MAD_survdec11","POR_survmay13")
survival_data <- pheno_data %>% 
  dplyr::select(site,block,pop,clon,tree,any_of(surv_measurements)) %>% 
  pivot_longer(cols=any_of(surv_measurements), names_to = NULL, values_to = "survival") %>% 
  drop_na(survival) 

survival_data[1:5,] %>% kable_mydf
site block pop clon tree survival
asturias 2 CAR CAR12 CAR12_2 1
asturias 4 STJ STJ14 STJ14_4 1
asturias 5 CUE CUE6 CUE6_5 1
asturias 6 MAD MAD1 MAD1_6 1
asturias 4 ORI ORI12 ORI12_4 1

3.2 Run the models

Stan code of the mortality model:

Code
stancode = stan_model(here("scripts/StanModels/ValidationCommonGarden_BinomialMortalityModel.stan"))
print(stancode)
S4 class stanmodel 'anon_model' coded as follows:
data {
  int N;
  int nb_tot[N];    // Total number of trees in the population
  int nb_dead[N];   // Number of dead trees in the population
  vector[N] H;      // Mean tree height of the population
  vector[N] X;      // Genomic offset or climatic transfer distance of the population
}
parameters {
  real beta_0;
  real beta_H;
  real beta_X;
}
model {
  nb_dead ~ binomial_logit(nb_tot,beta_0 + beta_H * H + beta_X * X);

  beta_0 ~ normal(0,5);//std_normal();
  beta_H ~ normal(0,5);//std_normal();
  beta_X ~ normal(0,5);//std_normal();
}
// generated quantities{
//   vector[N] log_lik;
//   // log likelihood for loo
//   for (n in 1:N) {
//     log_lik[n] = binomial_logit_lpmf( nb_dead[n] | nb_tot[n] , beta_0 + beta_H * H[n] + beta_X * X[n]);
//   }
// } 

We run the models for each common garden and each genomic offset predictions/climatic transfer distances, so a total of 5 * 15 = 75 models runs.

Code
coefftab <- lapply(unique(survival_data$site),function(site_i){
  
lapply(unique(df$method_input_code), function(method_input_code_i){

# Subset the data for the site i and method i
sub_data <- survival_data %>% 
  filter(site == site_i) %>% 
  group_by(pop) %>% 
  dplyr::summarise(nb_dead=n()-sum(survival),nb_tot=n()) # transform survival data into mortality data
    
sub_data <- df %>% 
  filter(method_input_code == method_input_code_i & cg == site_i) %>% 
  inner_join(sub_data, by="pop") %>% 
  inner_join(pop_heights %>% dplyr::select(any_of(c("height", "pop"))), by="pop") %>% 
  arrange(pop)
      
# Data in a list for Stan 
stanlist <- list(N = nrow(sub_data),
                 nb_dead = sub_data$nb_dead,
                 nb_tot = sub_data$nb_tot,
                 H = (sub_data$height-mean(sub_data$height)) / sd(sub_data$height),
                 X = (sub_data$varX -mean(sub_data$varX)) / sd(sub_data$varX))
    
# Running the model
mod <- sampling(stancode, data = stanlist, iter = 2000, chains = 4, cores = 4, save_warmup = FALSE) 
    
# Save the model and the stanlist
list(mod = mod, 
     stanlist = stanlist) %>% 
  saveRDS(file=here(paste0("outputs/ValidationCommonGarden/MortalityModels/stan_models/",
                                                           site_i,"_",method_input_code_i,".rds")))
    
  })
  
})
Code
coefftab <- lapply(unique(survival_data$site),function(site_i){
  
  lapply(unique(df$method_input_code), function(method_input_code_i){
    
    # Subset the data - keeping only one set of method / SNP set
    sub_data <- df %>% filter(method_input_code == method_input_code_i & cg == site_i)
  
    mod <- readRDS(file=here(paste0("outputs/ValidationCommonGarden/MortalityModels/stan_models/",
                                    site_i,"_",method_input_code_i,".rds")))[["mod"]]
                     
    # Save coefficients
    broom.mixed::tidyMCMC(mod,
                          droppars = NULL, 
                          robust = FALSE, # give mean and standard deviation
                          ess = F, 
                          rhat = F, 
                          conf.int = T,
                          conf.level = 0.95) %>% 
      filter(str_detect(term, c('beta'))) %>% 
      dplyr::rename(mean=estimate,
                    std_deviation=std.error,
                    conf_low=conf.low,
                    conf_high=conf.high) %>% 
      mutate(cg = site_i,
             method_input_code = method_input_code_i,
             method_input_name = unique(sub_data$method_input_name),
             method = unique(sub_data$method),
             input_name = unique(sub_data$input_name),
             input_code = unique(sub_data$input_code),
             .before=1)     

}) %>% bind_rows()
 
}) %>% bind_rows()
  
coefftab %>% saveRDS(file=here("outputs/ValidationCommonGarden/MortalityModels/coefftab.rds"))

3.3 Model coefficients

We load the model coefficients:

Code
coefftab <- readRDS(file=here("outputs/ValidationCommonGarden/MortalityModels/coefftab.rds"))
coefftab[1:5,] %>% kable_mydf()
cg method_input_code method_input_name method input_name input_code term mean std_deviation conf_low conf_high
asturias CTD_bio1 CTD - Mean annual temperature (bio1, °C) CTD Mean annual temperature (bio1, °C) bio1 beta_0 -2.99 0.07 -3.14 -2.85
asturias CTD_bio1 CTD - Mean annual temperature (bio1, °C) CTD Mean annual temperature (bio1, °C) bio1 beta_H 0.15 0.07 0.01 0.29
asturias CTD_bio1 CTD - Mean annual temperature (bio1, °C) CTD Mean annual temperature (bio1, °C) bio1 beta_X -0.02 0.08 -0.17 0.13
asturias CTD_bio12 CTD - Annual precipitation (bio12, mm) CTD Annual precipitation (bio12, mm) bio12 beta_0 -2.98 0.07 -3.12 -2.84
asturias CTD_bio12 CTD - Annual precipitation (bio12, mm) CTD Annual precipitation (bio12, mm) bio12 beta_H 0.13 0.08 -0.02 0.28

Below, we plot the mean and 95% credible intervals of:

  • the \(\beta_X\) coefficients, which stand for the association between the counts of dead trees and the genomic offset predictions (i.e. capturing the potential maladaptation of the populations when transplanted in the new environment of the common gardens) or the climatic transfer distances.

  • the \(\beta_H\) coefficients, which stand for the association between the counts of dead trees and the BLUPs for height in the five common gardens, used as a proxy of the initial seedling height (a confounder in the model).

Graph titles include the time in months corresponding to the age at which height and survival were recorded. Coefficients in the green area have the expected sign, reflecting higher mortality rates in populations with higher GO predictions.

Code
coeff_match <- list(beta_H = "Regression coefficients $\\beta_H$ in mortality models",
                    beta_X = "Regression coefficients $\\beta_X$ in mortality models")

p <- lapply(c("beta_X","beta_H"), function(coeff){
  
p <- coefftab %>% 
  filter(term == coeff) %>% 
  mutate(cg_name = case_when(cg == "asturias" ~ "Asturias, Spain (37 months)",
                             cg == "bordeaux" ~ "Bordeaux, France (85 months)",
                             cg == "caceres" ~ "Cáceres, Spain (8 months)",
                             cg == "madrid" ~ "Madrid, Spain (13 months)",
                             cg == "portugal" ~ "Fundão, Portugal (27 months)")) %>% 
  mutate(input_name = factor(input_name, levels = c(unique(coefftab$input_name)[[13]],
                                                    unique(coefftab$input_name)[[11]],
                                                    unique(coefftab$input_name)[[12]],
                                                    unique(coefftab$input_name)[[10]],
                                                    unique(coefftab$input_name)[[8]],
                                                    unique(coefftab$input_name)[[9]],
                                                    unique(coefftab$input_name)[[1]],
                                                    unique(coefftab$input_name)[[2]],
                                                    unique(coefftab$input_name)[[3]],
                                                    unique(coefftab$input_name)[[4]],
                                                    unique(coefftab$input_name)[[5]],
                                                    unique(coefftab$input_name)[[6]],
                                                    unique(coefftab$input_name)[[7]])),
         method = factor(method, levels = c(unique(coefftab$method)[[2]],
                                            unique(coefftab$method)[[3]],
                                            unique(coefftab$method)[[4]],
                                            unique(coefftab$method)[[5]],
                                            unique(coefftab$method)[[6]],
                                            unique(coefftab$method)[[1]])),
         cg_name = factor(cg_name, levels = c("Madrid, Spain (13 months)",
                                              "Cáceres, Spain (8 months)",
                                              "Asturias, Spain (37 months)",
                                              "Bordeaux, France (85 months)",
                                              "Fundão, Portugal (27 months)")))


# Plots with CTD and SNP sets
# ===========================
p_allvar <- p %>% ggplot(aes(x = method,
                             y = mean,
                             ymin = conf_low, 
                             ymax = conf_high,
                             color = input_name,
                             shape = input_name)) +
  {if(coeff=="beta_X")  geom_rect(inherit.aes=FALSE, aes(xmin=-Inf, xmax=Inf, ymin=0, ymax=Inf), color="transparent", fill="green", alpha=0.005)} + 
  geom_hline(yintercept = 0,color="gray") +
  geom_pointinterval(position = position_dodge(width = 0.6),
                     point_size=3, size=2) +
  facet_wrap(~cg_name, ncol=2) + 
  ylab(TeX(coeff_match[[coeff]])) +
  xlab("") +
  scale_color_manual(values=colors_coeff)+
  scale_shape_manual(values = c(rep(17,6),rep(16,7))) +
  theme_bw() +
  labs(color="",shape="") +
  theme(axis.text.x =  element_text(size=13),
        axis.text.y = element_text(size=13),
        axis.title.y = element_text(size=16),
        legend.title = element_text(size=13), 
        strip.text.x = element_text(size = 16),
        strip.background = element_blank(),
        legend.position = c(0.77,0.15),
        legend.text=element_text(size=12),
        panel.grid.minor.x=element_blank(),
        panel.grid.major.x=element_blank()) +
  guides(color = guide_legend(ncol=1),
         shape = guide_legend(override.aes = list(size =2 )))

# save in pdf and png
ggsave(p_allvar,
       file=here(paste0("figs/ValidationCommonGarden/MortalityModels/",coeff,"_SNPsandCTD.pdf")),
       device="pdf",
       height=11,
       width=13)

ggsave(p_allvar,
       file=here(paste0("figs/ValidationCommonGarden/MortalityModels/",coeff,"_SNPsandCTD.png")),
       height=11,
       width=13)


# Plots with SNPs only
# ====================
p_SNPs <- p %>% 
  filter(!method == "CTD") %>% 
  ggplot(aes(x = method,
             y = mean,
             ymin = conf_low, 
             ymax = conf_high,
             color = input_name,
             shape = input_name)) +
  {if(coeff=="beta_X")  geom_rect(inherit.aes=FALSE, aes(xmin=-Inf, xmax=Inf, ymin=0, ymax=Inf), color="transparent", fill="green", alpha=0.007)} + 
  geom_hline(yintercept = 0,color="gray") +
  geom_pointinterval(position = position_dodge(width = 0.45),
                     point_size=3, size=2) + # 
  facet_wrap(~cg_name, ncol=2) +
  ylab(TeX(coeff_match[[coeff]])) +
  xlab("") +
  scale_color_manual(values=colors_coeff) +
  scale_shape_manual(values = c(rep(17,6),rep(16,6))) +
  theme_bw() +
  labs(color="",shape="") +
  theme(axis.text.x =  element_text(size=13),
        axis.text.y = element_text(size=13),
        axis.title.y = element_text(size=16),
        axis.title.x = element_text(size=1),
        legend.title=element_text(size=13), 
        strip.text.x = element_text(size = 16),
        strip.background = element_blank(),
        legend.position = c(0.77,0.15),
        legend.text=element_text(size=14),
        panel.grid.minor.x=element_blank(),
        panel.grid.major.x=element_blank()) +
  guides(color=guide_legend(ncol=1),
         shape = guide_legend(override.aes = list(size =2 )))

# save in pdf and png
ggsave(p_SNPs,
       file=here(paste0("figs/ValidationCommonGarden/MortalityModels/",coeff,"_SNPsets.pdf")),
       device="pdf",
       height=11,
       width=13)

ggsave(p_SNPs,
       file=here(paste0("figs/ValidationCommonGarden/MortalityModels/",coeff,"_SNPsets.png")),
       height=11,
       width=13)

# Figure in the main manuscript
###############################

# We want to add some images to represent the climatic differences among common gardens

if(coeff=="beta_X"){
  
annotation_custom2 <- function (grob, xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = Inf, data){ layer(data = data, stat = StatIdentity, position = PositionIdentity, 
        geom = ggplot2:::GeomCustomAnn,
        inherit.aes = TRUE, params = list(grob = grob, 
                                          xmin = xmin, xmax = xmax, 
                                          ymin = ymin, ymax = ymax))}

df_png <- p %>% 
  dplyr::select(cg_name) %>%
  distinct %>% 
  dplyr::mutate(image = case_when(cg_name == "Asturias, Spain (37 months)" ~ "reports/cloud.png",
                                  cg_name == "Bordeaux, France (85 months)" ~ "reports/cloud.png",
                                  cg_name == "Cáceres, Spain (8 months)" ~ "reports/sun.png",
                                  cg_name == "Madrid, Spain (13 months)" ~ "reports/sun.png",
                                  cg_name == "Fundão, Portugal (27 months)" ~ 'reports/cloud_and_sun.png'))

list_pngs <- lapply(unique(p$cg_name), function(cg_name_i){
  
sub <-  p %>% 
  filter(cg_name==cg_name_i) %>%
  slice(1)

png_cg = annotation_custom2(rasterGrob(readPNG(here(df_png[df_png$cg_name==cg_name_i,"image"])),interpolate=TRUE), 
                            ymin = -0.54,
                            ymax= -0.24,
                            xmin = 0.42,
                            xmax = 1, 
                            data=sub)
})

p_allvar_images <- p_allvar + list_pngs[[1]]+ list_pngs[[2]]+ list_pngs[[3]]+ list_pngs[[4]]+ list_pngs[[5]]

# save in pdf and png
ggsave(p_allvar_images,
       file=here(paste0("figs/ValidationCommonGarden/MortalityModels/",coeff,"_SNPsets_WithCloudAndSunImages.pdf")),
       device="pdf",
       height=11,
       width=13)

}

return(p_allvar)

})
Warning: A numeric `legend.position` argument in `theme()` was deprecated in ggplot2 3.5.0.
ℹ Please use the `legend.position.inside` argument of `theme()` instead.
Warning in geom_rect(inherit.aes = FALSE, aes(xmin = -Inf, xmax = Inf, ymin = 0, : All aesthetics have length 1, but the data has 185 rows.
ℹ Please consider using `annotate()` or provide this layer with data containing a single row.
All aesthetics have length 1, but the data has 185 rows.
ℹ Please consider using `annotate()` or provide this layer with data containing a single row.
Warning in geom_rect(inherit.aes = FALSE, aes(xmin = -Inf, xmax = Inf, ymin = 0, : All aesthetics have length 1, but the data has 150 rows.
ℹ Please consider using `annotate()` or provide this layer with data containing a single row.
All aesthetics have length 1, but the data has 150 rows.
ℹ Please consider using `annotate()` or provide this layer with data containing a single row.
Warning in geom_rect(inherit.aes = FALSE, aes(xmin = -Inf, xmax = Inf, ymin = 0, : All aesthetics have length 1, but the data has 185 rows.
ℹ Please consider using `annotate()` or provide this layer with data containing a single row.
Code
p
[[1]]
Warning in geom_rect(inherit.aes = FALSE, aes(xmin = -Inf, xmax = Inf, ymin = 0, : All aesthetics have length 1, but the data has 185 rows.
ℹ Please consider using `annotate()` or provide this layer with data containing a single row.


[[2]]

Interpretation

As expected, the association between mortality rates and the initial tree height was particularly strong in Madrid and Cáceres (Spain) where the seedlings experienced an extreme drought event the same year they were planted, resulting in very high mortality rates. More surprisingly, this association was also osberved in Fundão, Portugal. However, the initial tree height was not associated with mortality rates in Bordeaux (France) and Asturias (Spain), which benefit from the favorable climates of the Atlantic region and in which the mortality rates were very low.

In Fundão, Cáceres and Madrid, for most genomic offset predictions, populations experiencing higher mortality rates also showed higher genomic offset. The most consistent associations across the three common gardens were obtained with:

  • Gradient Forest (GF) with both candidate and control SNPs.

  • Redundancy analysis (RDA) with both candidate and control SNPs.

  • Latent factor mixed model (LFMM) with all SNPs or control SNPs.

Interestingly, climatic transfer distances generally explain mortality probability less well than genomic offset predictions, and none of them show a consistent association with the number of dead trees across the three common gardens. Note that it is almost the case for the climatic transfer distance calculated based on the temperature seasonality, which shows a positive association with the counts of dead trees in Madrid and Fundão, and almost in Cáceres.

In Asturias, no association was detected between genomic offset predictions and mortality rates, which may be due to the favorable climatic conditions of this common garden and the associated low mortality rates (only 206 dead trees). Therefore, in this common garden, mortality events are likely to be mostly random events due to other factors than climate. However, we can note that an association was detected between mortality rates and isothermality: populations from areas with high isothermality tend to die more in the common garden of Asturias, Spain.

In Bordeaux, there were even less dead trees than in Asturias (119 dead trees). However, the two nonparametric approach used to predict the genomic offset (GDM and GF) detected a positive association between mortality rates and the genomic offset predictions obtained with both the candidate and the control SNPs. A positive association was also detected with the genomic offset predictions obtained with candidate SNPs with the LFMM approach. On the other hand, none of the genomic offset predictions from the RDA approach were associated with the mortality rates. Interestingly, the only association with a climatic transfer distance was with the temperature seasonality (bio4), which was the most important variable to explain the turnover in allele frequency in the GF and GDM approaches.

3.4 Predictions of tree mortality probability

3.4.1 Plots

Let’s visualize the uncertainty around the estimation of the probability of tree mortality \(p\).

Code
####################################################################################################
# Visualizing the relationships between GO predictions or CTDs and mortality probability predictions
####################################################################################################

coefftab <- readRDS(file = here("outputs/ValidationCommonGarden/MortalityModels/coefftab.rds")) %>%
  mutate(cg_name = case_when(cg == "asturias" ~ "Asturias, Spain (37 months)",
                             cg == "bordeaux" ~ "Bordeaux, France (85 months)",
                             cg == "caceres" ~ "Cáceres, Spain (8 months)",
                             cg == "madrid" ~ "Madrid, Spain (13 months)",
                             cg == "portugal" ~ "Fundão, Portugal (27 months)"))

lapply(unique(coefftab$cg), function(site_i){
  
  lapply(unique(coefftab$method_input_code), function(method_input_code_i){
    
    sub_coefftab <- coefftab %>% 
      filter(method_input_code == method_input_code_i & cg == site_i)
    
    # We load the stanlist
    stanlist <- readRDS(file = here(paste0("outputs/ValidationCommonGarden/MortalityModels/stan_models/",
                                       site_i,"_",method_input_code_i,".rds")))$stanlist

    # We load the model
    mod <- readRDS(file = here(paste0("outputs/ValidationCommonGarden/MortalityModels/stan_models/",
                                      site_i,"_",method_input_code_i,".rds")))$mod

    # We extract the posterior distributions of all parameters
    post <- as.data.frame(mod)

    # Vector of predictor values (based on the min and max of the GO predictions)
    x_vec <- seq(min(stanlist$X),max(stanlist$X),0.05)

    # Function to predict the mortality probability p with the initial tree height fixed to the mean
    f_p <- function(x) 1 / (exp(-(post$`beta_0` + post$`beta_H` * mean(stanlist$H) + post$`beta_X` * x)) + 1)

    p_pred <- 
      sapply(x_vec, f_p) %>% 
      tibble::as_tibble() %>%
      rename_all(function(x) x_vec) %>%
      mutate(Iter = row_number()) %>%
      gather(x, p, -Iter) %>%
      group_by(x) %>%
      mutate(hpdi_l = HDInterval::hdi(p, credMass = 0.90)[1],
             hpdi_h = HDInterval::hdi(p, credMass = 0.90)[2],
             p_mean = mean(p)) %>%
      ungroup() %>%
      mutate(x = as.numeric(x))
             
    # Keep mean and 90% HDPI of p (one value for each iteration)
    p_mean_df <- p_pred %>%
      dplyr::select(x,hpdi_l,hpdi_h,p_mean) %>% 
      dplyr::distinct()

    # Plots
    #=======

    # Plot options 
    y_limits <- c(0,1)
    if(unique(sub_coefftab$method) == "CTD"){x_axis <- "Mean-centered CTD"} else {x_axis <- "Mean-centered GO predictions"}
    y_axis <- "Probability of tree mortality"
    
    p <- ggplot() +
      ylim(y_limits) +
      labs(y=y_axis, x=x_axis) +
      theme_bw()

    # First 100 posterior draws of the mortality probability p
    p1 <- p +
      geom_point(data = p_pred %>% filter(Iter < 101),
                 aes(x, p), alpha = .2, color = 'dodgerblue')
  
   # Mean (line) and 90% HPDI (shaded region) of the mortality probability p
    p2 <- p +
      geom_ribbon(data = p_mean_df,
                  aes(x = x, ymin = hpdi_l, ymax = hpdi_h),
                  alpha = .2,
                  fill='dodgerblue') +
      geom_line(data = p_mean_df,
                aes(x=x, y=p_mean), col="dodgerblue4")

  p1_p2 <- ggarrange(p1, p2 + labs(y=""), nrow = 1)

  # Title
  title <- ggdraw() + 
    draw_label(paste0(unique(sub_coefftab$cg_name)," - ",unique(sub_coefftab$method_input_name)),
               fontface = 'bold',
               x = 0, 
               hjust = 0) +   
    theme(plot.margin = margin(0, 0, 0, 7))

  # merge title and plots 
  p1_p2 <- plot_grid(title, p1_p2, ncol = 1, rel_heights = c(0.1, 1))

ggsave(p1_p2,
       file=here(paste0("figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/",
                        method_input_code_i,"_",site_i,".pdf")),
       device="pdf",
       height=6,
       width=11)

})
})
[[1]]
[[1]][[1]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/CTD_bio1_asturias.pdf"

[[1]][[2]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/CTD_bio12_asturias.pdf"

[[1]][[3]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/CTD_bio15_asturias.pdf"

[[1]][[4]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/CTD_bio3_asturias.pdf"

[[1]][[5]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/CTD_bio4_asturias.pdf"

[[1]][[6]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/CTD_SHM_asturias.pdf"

[[1]][[7]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/CTD_EucliDist_asturias.pdf"

[[1]][[8]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/GDM_all_cand_asturias.pdf"

[[1]][[9]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/GF_all_cand_asturias.pdf"

[[1]][[10]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/LFMM_all_cand_asturias.pdf"

[[1]][[11]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/pRDA_all_cand_asturias.pdf"

[[1]][[12]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/RDA_all_cand_asturias.pdf"

[[1]][[13]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/GDM_common_cand_asturias.pdf"

[[1]][[14]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/GF_common_cand_asturias.pdf"

[[1]][[15]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/LFMM_common_cand_asturias.pdf"

[[1]][[16]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/pRDA_common_cand_asturias.pdf"

[[1]][[17]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/RDA_common_cand_asturias.pdf"

[[1]][[18]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/GDM_corrected_cand_asturias.pdf"

[[1]][[19]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/GF_corrected_cand_asturias.pdf"

[[1]][[20]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/LFMM_corrected_cand_asturias.pdf"

[[1]][[21]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/pRDA_corrected_cand_asturias.pdf"

[[1]][[22]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/RDA_corrected_cand_asturias.pdf"

[[1]][[23]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/GDM_randomfreq_control_asturias.pdf"

[[1]][[24]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/GF_randomfreq_control_asturias.pdf"

[[1]][[25]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/LFMM_randomfreq_control_asturias.pdf"

[[1]][[26]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/pRDA_randomfreq_control_asturias.pdf"

[[1]][[27]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/RDA_randomfreq_control_asturias.pdf"

[[1]][[28]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/GDM_simfreq_control_asturias.pdf"

[[1]][[29]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/GF_simfreq_control_asturias.pdf"

[[1]][[30]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/LFMM_simfreq_control_asturias.pdf"

[[1]][[31]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/pRDA_simfreq_control_asturias.pdf"

[[1]][[32]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/RDA_simfreq_control_asturias.pdf"

[[1]][[33]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/GDM_all_asturias.pdf"

[[1]][[34]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/GF_all_asturias.pdf"

[[1]][[35]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/LFMM_all_asturias.pdf"

[[1]][[36]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/pRDA_all_asturias.pdf"

[[1]][[37]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/RDA_all_asturias.pdf"


[[2]]
[[2]][[1]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/CTD_bio1_bordeaux.pdf"

[[2]][[2]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/CTD_bio12_bordeaux.pdf"

[[2]][[3]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/CTD_bio15_bordeaux.pdf"

[[2]][[4]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/CTD_bio3_bordeaux.pdf"

[[2]][[5]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/CTD_bio4_bordeaux.pdf"

[[2]][[6]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/CTD_SHM_bordeaux.pdf"

[[2]][[7]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/CTD_EucliDist_bordeaux.pdf"

[[2]][[8]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/GDM_all_cand_bordeaux.pdf"

[[2]][[9]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/GF_all_cand_bordeaux.pdf"

[[2]][[10]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/LFMM_all_cand_bordeaux.pdf"

[[2]][[11]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/pRDA_all_cand_bordeaux.pdf"

[[2]][[12]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/RDA_all_cand_bordeaux.pdf"

[[2]][[13]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/GDM_common_cand_bordeaux.pdf"

[[2]][[14]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/GF_common_cand_bordeaux.pdf"

[[2]][[15]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/LFMM_common_cand_bordeaux.pdf"

[[2]][[16]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/pRDA_common_cand_bordeaux.pdf"

[[2]][[17]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/RDA_common_cand_bordeaux.pdf"

[[2]][[18]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/GDM_corrected_cand_bordeaux.pdf"

[[2]][[19]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/GF_corrected_cand_bordeaux.pdf"

[[2]][[20]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/LFMM_corrected_cand_bordeaux.pdf"

[[2]][[21]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/pRDA_corrected_cand_bordeaux.pdf"

[[2]][[22]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/RDA_corrected_cand_bordeaux.pdf"

[[2]][[23]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/GDM_randomfreq_control_bordeaux.pdf"

[[2]][[24]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/GF_randomfreq_control_bordeaux.pdf"

[[2]][[25]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/LFMM_randomfreq_control_bordeaux.pdf"

[[2]][[26]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/pRDA_randomfreq_control_bordeaux.pdf"

[[2]][[27]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/RDA_randomfreq_control_bordeaux.pdf"

[[2]][[28]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/GDM_simfreq_control_bordeaux.pdf"

[[2]][[29]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/GF_simfreq_control_bordeaux.pdf"

[[2]][[30]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/LFMM_simfreq_control_bordeaux.pdf"

[[2]][[31]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/pRDA_simfreq_control_bordeaux.pdf"

[[2]][[32]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/RDA_simfreq_control_bordeaux.pdf"

[[2]][[33]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/GDM_all_bordeaux.pdf"

[[2]][[34]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/GF_all_bordeaux.pdf"

[[2]][[35]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/LFMM_all_bordeaux.pdf"

[[2]][[36]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/pRDA_all_bordeaux.pdf"

[[2]][[37]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/RDA_all_bordeaux.pdf"


[[3]]
[[3]][[1]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/CTD_bio1_caceres.pdf"

[[3]][[2]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/CTD_bio12_caceres.pdf"

[[3]][[3]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/CTD_bio15_caceres.pdf"

[[3]][[4]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/CTD_bio3_caceres.pdf"

[[3]][[5]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/CTD_bio4_caceres.pdf"

[[3]][[6]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/CTD_SHM_caceres.pdf"

[[3]][[7]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/CTD_EucliDist_caceres.pdf"

[[3]][[8]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/GDM_all_cand_caceres.pdf"

[[3]][[9]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/GF_all_cand_caceres.pdf"

[[3]][[10]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/LFMM_all_cand_caceres.pdf"

[[3]][[11]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/pRDA_all_cand_caceres.pdf"

[[3]][[12]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/RDA_all_cand_caceres.pdf"

[[3]][[13]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/GDM_common_cand_caceres.pdf"

[[3]][[14]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/GF_common_cand_caceres.pdf"

[[3]][[15]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/LFMM_common_cand_caceres.pdf"

[[3]][[16]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/pRDA_common_cand_caceres.pdf"

[[3]][[17]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/RDA_common_cand_caceres.pdf"

[[3]][[18]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/GDM_corrected_cand_caceres.pdf"

[[3]][[19]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/GF_corrected_cand_caceres.pdf"

[[3]][[20]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/LFMM_corrected_cand_caceres.pdf"

[[3]][[21]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/pRDA_corrected_cand_caceres.pdf"

[[3]][[22]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/RDA_corrected_cand_caceres.pdf"

[[3]][[23]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/GDM_randomfreq_control_caceres.pdf"

[[3]][[24]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/GF_randomfreq_control_caceres.pdf"

[[3]][[25]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/LFMM_randomfreq_control_caceres.pdf"

[[3]][[26]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/pRDA_randomfreq_control_caceres.pdf"

[[3]][[27]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/RDA_randomfreq_control_caceres.pdf"

[[3]][[28]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/GDM_simfreq_control_caceres.pdf"

[[3]][[29]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/GF_simfreq_control_caceres.pdf"

[[3]][[30]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/LFMM_simfreq_control_caceres.pdf"

[[3]][[31]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/pRDA_simfreq_control_caceres.pdf"

[[3]][[32]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/RDA_simfreq_control_caceres.pdf"

[[3]][[33]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/GDM_all_caceres.pdf"

[[3]][[34]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/GF_all_caceres.pdf"

[[3]][[35]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/LFMM_all_caceres.pdf"

[[3]][[36]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/pRDA_all_caceres.pdf"

[[3]][[37]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/RDA_all_caceres.pdf"


[[4]]
[[4]][[1]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/CTD_bio1_madrid.pdf"

[[4]][[2]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/CTD_bio12_madrid.pdf"

[[4]][[3]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/CTD_bio15_madrid.pdf"

[[4]][[4]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/CTD_bio3_madrid.pdf"

[[4]][[5]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/CTD_bio4_madrid.pdf"

[[4]][[6]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/CTD_SHM_madrid.pdf"

[[4]][[7]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/CTD_EucliDist_madrid.pdf"

[[4]][[8]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/GDM_all_cand_madrid.pdf"

[[4]][[9]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/GF_all_cand_madrid.pdf"

[[4]][[10]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/LFMM_all_cand_madrid.pdf"

[[4]][[11]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/pRDA_all_cand_madrid.pdf"

[[4]][[12]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/RDA_all_cand_madrid.pdf"

[[4]][[13]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/GDM_common_cand_madrid.pdf"

[[4]][[14]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/GF_common_cand_madrid.pdf"

[[4]][[15]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/LFMM_common_cand_madrid.pdf"

[[4]][[16]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/pRDA_common_cand_madrid.pdf"

[[4]][[17]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/RDA_common_cand_madrid.pdf"

[[4]][[18]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/GDM_corrected_cand_madrid.pdf"

[[4]][[19]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/GF_corrected_cand_madrid.pdf"

[[4]][[20]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/LFMM_corrected_cand_madrid.pdf"

[[4]][[21]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/pRDA_corrected_cand_madrid.pdf"

[[4]][[22]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/RDA_corrected_cand_madrid.pdf"

[[4]][[23]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/GDM_randomfreq_control_madrid.pdf"

[[4]][[24]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/GF_randomfreq_control_madrid.pdf"

[[4]][[25]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/LFMM_randomfreq_control_madrid.pdf"

[[4]][[26]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/pRDA_randomfreq_control_madrid.pdf"

[[4]][[27]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/RDA_randomfreq_control_madrid.pdf"

[[4]][[28]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/GDM_simfreq_control_madrid.pdf"

[[4]][[29]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/GF_simfreq_control_madrid.pdf"

[[4]][[30]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/LFMM_simfreq_control_madrid.pdf"

[[4]][[31]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/pRDA_simfreq_control_madrid.pdf"

[[4]][[32]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/RDA_simfreq_control_madrid.pdf"

[[4]][[33]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/GDM_all_madrid.pdf"

[[4]][[34]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/GF_all_madrid.pdf"

[[4]][[35]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/LFMM_all_madrid.pdf"

[[4]][[36]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/pRDA_all_madrid.pdf"

[[4]][[37]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/RDA_all_madrid.pdf"


[[5]]
[[5]][[1]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/CTD_bio1_portugal.pdf"

[[5]][[2]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/CTD_bio12_portugal.pdf"

[[5]][[3]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/CTD_bio15_portugal.pdf"

[[5]][[4]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/CTD_bio3_portugal.pdf"

[[5]][[5]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/CTD_bio4_portugal.pdf"

[[5]][[6]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/CTD_SHM_portugal.pdf"

[[5]][[7]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/CTD_EucliDist_portugal.pdf"

[[5]][[8]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/GDM_all_cand_portugal.pdf"

[[5]][[9]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/GF_all_cand_portugal.pdf"

[[5]][[10]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/LFMM_all_cand_portugal.pdf"

[[5]][[11]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/pRDA_all_cand_portugal.pdf"

[[5]][[12]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/RDA_all_cand_portugal.pdf"

[[5]][[13]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/GDM_common_cand_portugal.pdf"

[[5]][[14]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/GF_common_cand_portugal.pdf"

[[5]][[15]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/LFMM_common_cand_portugal.pdf"

[[5]][[16]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/pRDA_common_cand_portugal.pdf"

[[5]][[17]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/RDA_common_cand_portugal.pdf"

[[5]][[18]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/GDM_corrected_cand_portugal.pdf"

[[5]][[19]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/GF_corrected_cand_portugal.pdf"

[[5]][[20]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/LFMM_corrected_cand_portugal.pdf"

[[5]][[21]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/pRDA_corrected_cand_portugal.pdf"

[[5]][[22]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/RDA_corrected_cand_portugal.pdf"

[[5]][[23]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/GDM_randomfreq_control_portugal.pdf"

[[5]][[24]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/GF_randomfreq_control_portugal.pdf"

[[5]][[25]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/LFMM_randomfreq_control_portugal.pdf"

[[5]][[26]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/pRDA_randomfreq_control_portugal.pdf"

[[5]][[27]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/RDA_randomfreq_control_portugal.pdf"

[[5]][[28]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/GDM_simfreq_control_portugal.pdf"

[[5]][[29]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/GF_simfreq_control_portugal.pdf"

[[5]][[30]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/LFMM_simfreq_control_portugal.pdf"

[[5]][[31]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/pRDA_simfreq_control_portugal.pdf"

[[5]][[32]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/RDA_simfreq_control_portugal.pdf"

[[5]][[33]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/GDM_all_portugal.pdf"

[[5]][[34]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/GF_all_portugal.pdf"

[[5]][[35]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/LFMM_all_portugal.pdf"

[[5]][[36]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/pRDA_all_portugal.pdf"

[[5]][[37]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/MortalityProbabilityPredictions_UncertaintyIntervals/RDA_all_portugal.pdf"

3.4.2 Table

In the table below:

  • mean and 90% credible intervals (highest posterior density intervals) of \(p\) for \(x=\{-1, 0, 1\}\).

  • percent change in \(p\) associated with a one-unit increase in \(x\) (which corresponds to a one-standard deviation increase).

Code
p_pred <- lapply(names(cg_names), function(site_i){
  
lapply(unique(df$method_input_code), function(method_input_code_i){   
#lapply(names(methods), function(method_i){

# We load the stanlist
stanlist <- readRDS(file = here(paste0("outputs/ValidationCommonGarden/MortalityModels/stan_models/",
                                       site_i,"_",method_input_code_i,".rds")))$stanlist

# We load the model
mod <- readRDS(file = here(paste0("outputs/ValidationCommonGarden/MortalityModels/stan_models/",
                                  site_i,"_",method_input_code_i,".rds")))$mod

# We extract the posterior distributions of all parameters
post <- as.data.frame(mod)

# Vector of predictor values (based on the min and max of the GO predictions)
x_vec <- c(-1,0,1)

# Function to predict the mortality probability p with the initial tree height fixed to the mean
f_p <- function(x) 1 / (exp(-(post$`beta_0` + post$`beta_H` * mean(stanlist$H) + post$`beta_X` * x)) + 1)

p_pred <- 
  sapply(x_vec, f_p) %>% 
  tibble::as_tibble() %>%
  rename_all(function(x) x_vec) %>%
  mutate(Iter = row_number()) %>%
  gather(x, p, -Iter) %>%
  group_by(x) %>%
  mutate(hpdi_l = HDInterval::hdi(p, credMass = 0.90)[1], # Highest posterior density intervals
         hpdi_h = HDInterval::hdi(p, credMass = 0.90)[2]
         # pi_l = rethinking::PI(p, prob = 0.90)[1], # Highest posterior percentile intervals
         # pi_h = rethinking::PI(p, prob = 0.90)[2]
         ) %>%
  mutate(p_mean = mean(p)) %>%
  ungroup() %>%
  mutate(x = as.numeric(x)) %>% 
  dplyr::select(x,p_mean,hpdi_l,hpdi_h) %>% 
  distinct() %>% 
  round(3)

delta_p_pos <- (p_pred$p_mean[p_pred$x==1] * 100 / p_pred$p_mean[p_pred$x==0])-100
delta_p_neg <- (p_pred$p_mean[p_pred$x==-1] * 100 / p_pred$p_mean[p_pred$x==0])-100

tibble(Site = cg_names[[site_i]],
       "Method" = df %>% filter(method_input_code == method_input_code_i) %>% pull(method_input_name) %>% unique(),
       "p(x=-1)"=paste0(p_pred$p_mean[p_pred$x==-1], " [",
                       p_pred$hpdi_l[p_pred$x==-1],";",
                       p_pred$hpdi_h[p_pred$x==-1],"]"),
       "p(x=0)"=paste0(p_pred$p_mean[p_pred$x==0], " [",
                       p_pred$hpdi_l[p_pred$x==0],";",
                       p_pred$hpdi_h[p_pred$x==0],"]"),
       "p(x=1)"=paste0(p_pred$p_mean[p_pred$x==1], " [",
                       p_pred$hpdi_l[p_pred$x==1],";",
                       p_pred$hpdi_h[p_pred$x==1],"]"),
       "+Δp"=paste0(round(delta_p_pos,3),"%"),
       "-Δp"=paste0(round(delta_p_neg,3),"%"))

}) %>% bind_rows()

}) %>% bind_rows()


p_pred %>% saveRDS(file = here("tables/MortalityPredCGs.rds"))
Code
readRDS(file = here("tables/MortalityPredCGs.rds"))  %>% kable_mydf()
Site Method p(x=-1) p(x=0) p(x=1) +Δp -Δp
Cáceres CTD - Mean annual temperature (bio1, °C) 0.915 [0.905;0.925] 0.916 [0.909;0.923] 0.917 [0.906;0.929] 0.109% -0.109%
Cáceres CTD - Annual precipitation (bio12, mm) 0.911 [0.9;0.921] 0.917 [0.91;0.924] 0.922 [0.911;0.933] 0.545% -0.654%
Cáceres CTD - Precipitation seasonality (bio15, index) 0.909 [0.897;0.922] 0.915 [0.908;0.922] 0.92 [0.911;0.931] 0.546% -0.656%
Cáceres CTD - Isothermality (bio3, index) 0.919 [0.908;0.93] 0.916 [0.908;0.923] 0.913 [0.902;0.923] -0.328% 0.328%
Cáceres CTD - Temperature seasonality (bio4, °C) 0.905 [0.892;0.917] 0.917 [0.91;0.924] 0.927 [0.916;0.938] 1.091% -1.309%
Cáceres CTD - Summer heat moisture index (SHM, °C/mm) 0.904 [0.893;0.917] 0.916 [0.909;0.923] 0.927 [0.916;0.938] 1.201% -1.31%
Cáceres CTD - All climatic variables (Euclidean distance) 0.908 [0.897;0.918] 0.916 [0.909;0.924] 0.924 [0.914;0.934] 0.873% -0.873%
Cáceres GDM - All candidate SNPs (380) 0.905 [0.893;0.919] 0.916 [0.909;0.923] 0.926 [0.914;0.937] 1.092% -1.201%
Cáceres GF - All candidate SNPs (380) 0.902 [0.89;0.915] 0.917 [0.91;0.924] 0.929 [0.919;0.94] 1.309% -1.636%
Cáceres LFMM - All candidate SNPs (380) 0.902 [0.887;0.915] 0.916 [0.909;0.923] 0.928 [0.917;0.938] 1.31% -1.528%
Cáceres pRDA - All candidate SNPs (380) 0.9 [0.887;0.914] 0.917 [0.91;0.924] 0.931 [0.92;0.942] 1.527% -1.854%
Cáceres RDA - All candidate SNPs (380) 0.898 [0.885;0.913] 0.916 [0.91;0.924] 0.931 [0.92;0.942] 1.638% -1.965%
Cáceres GDM - Common candidate SNPs (69) 0.906 [0.895;0.918] 0.917 [0.91;0.925] 0.926 [0.916;0.938] 0.981% -1.2%
Cáceres GF - Common candidate SNPs (69) 0.904 [0.892;0.915] 0.917 [0.91;0.924] 0.928 [0.917;0.939] 1.2% -1.418%
Cáceres LFMM - Common candidate SNPs (69) 0.908 [0.895;0.921] 0.916 [0.908;0.923] 0.922 [0.912;0.933] 0.655% -0.873%
Cáceres pRDA - Common candidate SNPs (69) 0.926 [0.916;0.936] 0.916 [0.909;0.923] 0.905 [0.893;0.916] -1.201% 1.092%
Cáceres RDA - Common candidate SNPs (69) 0.899 [0.887;0.913] 0.918 [0.909;0.924] 0.932 [0.921;0.943] 1.525% -2.07%
Cáceres GDM - Candidate SNPs considering pop. struct. correction (221) 0.902 [0.887;0.916] 0.916 [0.909;0.923] 0.928 [0.917;0.94] 1.31% -1.528%
Cáceres GF - Candidate SNPs considering pop. struct. correction (221) 0.901 [0.887;0.914] 0.917 [0.909;0.924] 0.93 [0.919;0.941] 1.418% -1.745%
Cáceres LFMM - Candidate SNPs considering pop. struct. correction (221) 0.904 [0.891;0.917] 0.916 [0.908;0.923] 0.926 [0.915;0.937] 1.092% -1.31%
Cáceres pRDA - Candidate SNPs considering pop. struct. correction (221) 0.901 [0.888;0.915] 0.917 [0.909;0.924] 0.929 [0.918;0.94] 1.309% -1.745%
Cáceres RDA - Candidate SNPs considering pop. struct. correction (221) 0.897 [0.881;0.911] 0.916 [0.909;0.923] 0.932 [0.921;0.943] 1.747% -2.074%
Cáceres GDM - Control SNPs unmatching allele frequencies (380) 0.909 [0.893;0.926] 0.915 [0.908;0.923] 0.92 [0.908;0.933] 0.546% -0.656%
Cáceres GF - Control SNPs unmatching allele frequencies (380) 0.9 [0.886;0.915] 0.915 [0.908;0.923] 0.928 [0.917;0.939] 1.421% -1.639%
Cáceres LFMM - Control SNPs unmatching allele frequencies (380) 0.903 [0.89;0.916] 0.916 [0.909;0.923] 0.927 [0.916;0.938] 1.201% -1.419%
Cáceres pRDA - Control SNPs unmatching allele frequencies (380) 0.9 [0.886;0.913] 0.917 [0.909;0.923] 0.931 [0.919;0.942] 1.527% -1.854%
Cáceres RDA - Control SNPs unmatching allele frequencies (380) 0.903 [0.889;0.918] 0.915 [0.907;0.922] 0.925 [0.914;0.935] 1.093% -1.311%
Cáceres GDM - Control SNPs matching allele frequencies (380) 0.909 [0.895;0.923] 0.916 [0.908;0.922] 0.921 [0.91;0.934] 0.546% -0.764%
Cáceres GF - Control SNPs matching allele frequencies (380) 0.9 [0.886;0.915] 0.916 [0.909;0.923] 0.929 [0.917;0.939] 1.419% -1.747%
Cáceres LFMM - Control SNPs matching allele frequencies (380) 0.903 [0.889;0.916] 0.916 [0.908;0.923] 0.927 [0.916;0.938] 1.201% -1.419%
Cáceres pRDA - Control SNPs matching allele frequencies (380) 0.897 [0.883;0.912] 0.916 [0.909;0.923] 0.932 [0.921;0.943] 1.747% -2.074%
Cáceres RDA - Control SNPs matching allele frequencies (380) 0.901 [0.885;0.915] 0.915 [0.908;0.922] 0.927 [0.917;0.937] 1.311% -1.53%
Cáceres GDM - All SNPs (9817) 0.906 [0.89;0.921] 0.915 [0.908;0.922] 0.923 [0.911;0.935] 0.874% -0.984%
Cáceres GF - All SNPs (9817) 0.9 [0.885;0.915] 0.916 [0.909;0.923] 0.929 [0.918;0.94] 1.419% -1.747%
Cáceres LFMM - All SNPs (9817) 0.902 [0.889;0.916] 0.916 [0.909;0.923] 0.928 [0.917;0.94] 1.31% -1.528%
Cáceres pRDA - All SNPs (9817) 0.899 [0.885;0.914] 0.917 [0.91;0.924] 0.931 [0.921;0.943] 1.527% -1.963%
Cáceres RDA - All SNPs (9817) 0.9 [0.885;0.916] 0.915 [0.908;0.923] 0.927 [0.917;0.938] 1.311% -1.639%
Asturias CTD - Mean annual temperature (bio1, °C) 0.049 [0.041;0.057] 0.048 [0.043;0.054] 0.047 [0.04;0.056] -2.083% 2.083%
Asturias CTD - Annual precipitation (bio12, mm) 0.051 [0.041;0.06] 0.048 [0.043;0.054] 0.047 [0.039;0.054] -2.083% 6.25%
Asturias CTD - Precipitation seasonality (bio15, index) 0.043 [0.035;0.051] 0.048 [0.043;0.054] 0.055 [0.045;0.065] 14.583% -10.417%
Asturias CTD - Isothermality (bio3, index) 0.038 [0.031;0.046] 0.047 [0.042;0.053] 0.058 [0.049;0.067] 23.404% -19.149%
Asturias CTD - Temperature seasonality (bio4, °C) 0.058 [0.046;0.069] 0.048 [0.043;0.054] 0.04 [0.031;0.048] -16.667% 20.833%
Asturias CTD - Summer heat moisture index (SHM, °C/mm) 0.053 [0.044;0.063] 0.048 [0.042;0.053] 0.043 [0.034;0.052] -10.417% 10.417%
Asturias CTD - All climatic variables (Euclidean distance) 0.044 [0.036;0.052] 0.048 [0.042;0.053] 0.052 [0.044;0.06] 8.333% -8.333%
Asturias GDM - All candidate SNPs (380) 0.057 [0.048;0.067] 0.047 [0.042;0.053] 0.04 [0.031;0.048] -14.894% 21.277%
Asturias GF - All candidate SNPs (380) 0.056 [0.046;0.067] 0.048 [0.043;0.054] 0.042 [0.034;0.05] -12.5% 16.667%
Asturias LFMM - All candidate SNPs (380) 0.053 [0.042;0.062] 0.048 [0.043;0.054] 0.044 [0.035;0.053] -8.333% 10.417%
Asturias pRDA - All candidate SNPs (380) 0.059 [0.048;0.069] 0.047 [0.042;0.053] 0.039 [0.03;0.047] -17.021% 25.532%
Asturias RDA - All candidate SNPs (380) 0.05 [0.041;0.06] 0.048 [0.043;0.054] 0.047 [0.038;0.055] -2.083% 4.167%
Asturias GDM - Common candidate SNPs (69) 0.059 [0.048;0.068] 0.048 [0.042;0.053] 0.039 [0.031;0.047] -18.75% 22.917%
Asturias GF - Common candidate SNPs (69) 0.057 [0.048;0.068] 0.048 [0.043;0.054] 0.041 [0.033;0.048] -14.583% 18.75%
Asturias LFMM - Common candidate SNPs (69) 0.052 [0.042;0.061] 0.048 [0.043;0.054] 0.045 [0.036;0.053] -6.25% 8.333%
Asturias pRDA - Common candidate SNPs (69) 0.052 [0.043;0.062] 0.048 [0.043;0.054] 0.045 [0.036;0.054] -6.25% 8.333%
Asturias RDA - Common candidate SNPs (69) 0.053 [0.043;0.063] 0.048 [0.042;0.053] 0.044 [0.035;0.053] -8.333% 10.417%
Asturias GDM - Candidate SNPs considering pop. struct. correction (221) 0.059 [0.047;0.069] 0.047 [0.042;0.053] 0.039 [0.029;0.047] -17.021% 25.532%
Asturias GF - Candidate SNPs considering pop. struct. correction (221) 0.057 [0.046;0.068] 0.048 [0.042;0.053] 0.041 [0.032;0.049] -14.583% 18.75%
Asturias LFMM - Candidate SNPs considering pop. struct. correction (221) 0.054 [0.043;0.064] 0.048 [0.042;0.054] 0.043 [0.033;0.052] -10.417% 12.5%
Asturias pRDA - Candidate SNPs considering pop. struct. correction (221) 0.058 [0.047;0.068] 0.047 [0.042;0.053] 0.039 [0.031;0.048] -17.021% 23.404%
Asturias RDA - Candidate SNPs considering pop. struct. correction (221) 0.056 [0.045;0.066] 0.048 [0.042;0.053] 0.041 [0.032;0.05] -14.583% 16.667%
Asturias GDM - Control SNPs unmatching allele frequencies (380) 0.046 [0.037;0.054] 0.049 [0.043;0.054] 0.052 [0.04;0.064] 6.122% -6.122%
Asturias GF - Control SNPs unmatching allele frequencies (380) 0.054 [0.043;0.064] 0.048 [0.042;0.054] 0.043 [0.034;0.052] -10.417% 12.5%
Asturias LFMM - Control SNPs unmatching allele frequencies (380) 0.047 [0.039;0.056] 0.048 [0.043;0.054] 0.049 [0.04;0.057] 2.083% -2.083%
Asturias pRDA - Control SNPs unmatching allele frequencies (380) 0.057 [0.046;0.067] 0.047 [0.042;0.053] 0.04 [0.031;0.05] -14.894% 21.277%
Asturias RDA - Control SNPs unmatching allele frequencies (380) 0.044 [0.036;0.052] 0.048 [0.042;0.053] 0.052 [0.043;0.06] 8.333% -8.333%
Asturias GDM - Control SNPs matching allele frequencies (380) 0.05 [0.04;0.058] 0.048 [0.042;0.053] 0.046 [0.036;0.057] -4.167% 4.167%
Asturias GF - Control SNPs matching allele frequencies (380) 0.054 [0.043;0.065] 0.048 [0.043;0.054] 0.043 [0.034;0.052] -10.417% 12.5%
Asturias LFMM - Control SNPs matching allele frequencies (380) 0.046 [0.038;0.054] 0.048 [0.042;0.054] 0.05 [0.042;0.059] 4.167% -4.167%
Asturias pRDA - Control SNPs matching allele frequencies (380) 0.056 [0.047;0.066] 0.047 [0.042;0.053] 0.04 [0.031;0.048] -14.894% 19.149%
Asturias RDA - Control SNPs matching allele frequencies (380) 0.046 [0.038;0.054] 0.048 [0.043;0.054] 0.051 [0.043;0.06] 6.25% -4.167%
Asturias GDM - All SNPs (9817) 0.05 [0.04;0.059] 0.048 [0.042;0.054] 0.046 [0.034;0.058] -4.167% 4.167%
Asturias GF - All SNPs (9817) 0.054 [0.044;0.065] 0.048 [0.043;0.054] 0.043 [0.034;0.051] -10.417% 12.5%
Asturias LFMM - All SNPs (9817) 0.048 [0.039;0.056] 0.048 [0.043;0.054] 0.049 [0.04;0.057] 2.083% 0%
Asturias pRDA - All SNPs (9817) 0.06 [0.049;0.071] 0.047 [0.041;0.052] 0.037 [0.028;0.046] -21.277% 27.66%
Asturias RDA - All SNPs (9817) 0.046 [0.037;0.054] 0.048 [0.043;0.054] 0.051 [0.043;0.059] 6.25% -4.167%
Madrid CTD - Mean annual temperature (bio1, °C) 0.761 [0.746;0.776] 0.752 [0.741;0.763] 0.743 [0.724;0.761] -1.197% 1.197%
Madrid CTD - Annual precipitation (bio12, mm) 0.738 [0.723;0.754] 0.757 [0.746;0.768] 0.774 [0.757;0.791] 2.246% -2.51%
Madrid CTD - Precipitation seasonality (bio15, index) 0.743 [0.728;0.76] 0.754 [0.743;0.765] 0.763 [0.747;0.779] 1.194% -1.459%
Madrid CTD - Isothermality (bio3, index) 0.731 [0.711;0.75] 0.752 [0.74;0.762] 0.772 [0.757;0.787] 2.66% -2.793%
Madrid CTD - Temperature seasonality (bio4, °C) 0.712 [0.694;0.733] 0.757 [0.747;0.769] 0.797 [0.779;0.814] 5.284% -5.945%
Madrid CTD - Summer heat moisture index (SHM, °C/mm) 0.742 [0.724;0.76] 0.754 [0.743;0.765] 0.766 [0.749;0.784] 1.592% -1.592%
Madrid CTD - All climatic variables (Euclidean distance) 0.733 [0.716;0.748] 0.757 [0.746;0.768] 0.78 [0.761;0.796] 3.038% -3.17%
Madrid GDM - All candidate SNPs (380) 0.721 [0.704;0.739] 0.758 [0.747;0.769] 0.792 [0.776;0.81] 4.485% -4.881%
Madrid GF - All candidate SNPs (380) 0.72 [0.701;0.738] 0.758 [0.747;0.769] 0.792 [0.774;0.809] 4.485% -5.013%
Madrid LFMM - All candidate SNPs (380) 0.725 [0.707;0.745] 0.756 [0.745;0.767] 0.783 [0.766;0.803] 3.571% -4.101%
Madrid pRDA - All candidate SNPs (380) 0.715 [0.695;0.733] 0.756 [0.745;0.767] 0.794 [0.776;0.812] 5.026% -5.423%
Madrid RDA - All candidate SNPs (380) 0.713 [0.694;0.731] 0.758 [0.747;0.768] 0.797 [0.778;0.813] 5.145% -5.937%
Madrid GDM - Common candidate SNPs (69) 0.718 [0.7;0.735] 0.758 [0.747;0.768] 0.794 [0.776;0.81] 4.749% -5.277%
Madrid GF - Common candidate SNPs (69) 0.72 [0.702;0.737] 0.758 [0.748;0.77] 0.792 [0.775;0.809] 4.485% -5.013%
Madrid LFMM - Common candidate SNPs (69) 0.734 [0.715;0.754] 0.754 [0.743;0.765] 0.773 [0.754;0.79] 2.52% -2.653%
Madrid pRDA - Common candidate SNPs (69) 0.745 [0.729;0.759] 0.755 [0.743;0.766] 0.765 [0.748;0.782] 1.325% -1.325%
Madrid RDA - Common candidate SNPs (69) 0.715 [0.698;0.735] 0.758 [0.746;0.769] 0.796 [0.78;0.814] 5.013% -5.673%
Madrid GDM - Candidate SNPs considering pop. struct. correction (221) 0.721 [0.702;0.738] 0.758 [0.746;0.769] 0.79 [0.773;0.809] 4.222% -4.881%
Madrid GF - Candidate SNPs considering pop. struct. correction (221) 0.717 [0.699;0.737] 0.757 [0.746;0.768] 0.792 [0.774;0.81] 4.624% -5.284%
Madrid LFMM - Candidate SNPs considering pop. struct. correction (221) 0.732 [0.713;0.752] 0.754 [0.743;0.765] 0.776 [0.757;0.794] 2.918% -2.918%
Madrid pRDA - Candidate SNPs considering pop. struct. correction (221) 0.719 [0.698;0.738] 0.756 [0.745;0.767] 0.79 [0.771;0.806] 4.497% -4.894%
Madrid RDA - Candidate SNPs considering pop. struct. correction (221) 0.707 [0.688;0.726] 0.756 [0.746;0.767] 0.8 [0.784;0.819] 5.82% -6.481%
Madrid GDM - Control SNPs unmatching allele frequencies (380) 0.735 [0.717;0.754] 0.759 [0.747;0.771] 0.781 [0.759;0.808] 2.899% -3.162%
Madrid GF - Control SNPs unmatching allele frequencies (380) 0.723 [0.705;0.741] 0.758 [0.747;0.769] 0.789 [0.771;0.807] 4.09% -4.617%
Madrid LFMM - Control SNPs unmatching allele frequencies (380) 0.734 [0.716;0.752] 0.756 [0.745;0.767] 0.776 [0.758;0.795] 2.646% -2.91%
Madrid pRDA - Control SNPs unmatching allele frequencies (380) 0.716 [0.697;0.737] 0.756 [0.744;0.766] 0.792 [0.774;0.81] 4.762% -5.291%
Madrid RDA - Control SNPs unmatching allele frequencies (380) 0.711 [0.693;0.73] 0.755 [0.743;0.765] 0.793 [0.777;0.811] 5.033% -5.828%
Madrid GDM - Control SNPs matching allele frequencies (380) 0.734 [0.719;0.749] 0.76 [0.749;0.771] 0.783 [0.766;0.802] 3.026% -3.421%
Madrid GF - Control SNPs matching allele frequencies (380) 0.722 [0.704;0.74] 0.758 [0.747;0.769] 0.789 [0.771;0.807] 4.09% -4.749%
Madrid LFMM - Control SNPs matching allele frequencies (380) 0.734 [0.715;0.753] 0.755 [0.744;0.766] 0.774 [0.755;0.793] 2.517% -2.781%
Madrid pRDA - Control SNPs matching allele frequencies (380) 0.711 [0.689;0.733] 0.755 [0.743;0.766] 0.793 [0.775;0.811] 5.033% -5.828%
Madrid RDA - Control SNPs matching allele frequencies (380) 0.711 [0.691;0.729] 0.754 [0.743;0.764] 0.792 [0.776;0.808] 5.04% -5.703%
Madrid GDM - All SNPs (9817) 0.734 [0.718;0.75] 0.76 [0.748;0.771] 0.784 [0.765;0.802] 3.158% -3.421%
Madrid GF - All SNPs (9817) 0.722 [0.704;0.74] 0.758 [0.746;0.769] 0.79 [0.773;0.81] 4.222% -4.749%
Madrid LFMM - All SNPs (9817) 0.73 [0.712;0.75] 0.755 [0.745;0.767] 0.779 [0.761;0.798] 3.179% -3.311%
Madrid pRDA - All SNPs (9817) 0.714 [0.695;0.733] 0.756 [0.745;0.767] 0.793 [0.775;0.811] 4.894% -5.556%
Madrid RDA - All SNPs (9817) 0.71 [0.691;0.729] 0.755 [0.744;0.767] 0.795 [0.779;0.812] 5.298% -5.96%
Fundão CTD - Mean annual temperature (bio1, °C) 0.39 [0.37;0.408] 0.361 [0.35;0.374] 0.334 [0.315;0.352] -7.479% 8.033%
Fundão CTD - Annual precipitation (bio12, mm) 0.4 [0.379;0.423] 0.369 [0.356;0.381] 0.339 [0.323;0.355] -8.13% 8.401%
Fundão CTD - Precipitation seasonality (bio15, index) 0.347 [0.328;0.368] 0.361 [0.349;0.374] 0.375 [0.356;0.391] 3.878% -3.878%
Fundão CTD - Isothermality (bio3, index) 0.346 [0.325;0.365] 0.363 [0.35;0.374] 0.38 [0.361;0.399] 4.683% -4.683%
Fundão CTD - Temperature seasonality (bio4, °C) 0.325 [0.305;0.343] 0.366 [0.353;0.378] 0.41 [0.39;0.433] 12.022% -11.202%
Fundão CTD - Summer heat moisture index (SHM, °C/mm) 0.325 [0.308;0.344] 0.366 [0.353;0.378] 0.409 [0.39;0.429] 11.749% -11.202%
Fundão CTD - All climatic variables (Euclidean distance) 0.34 [0.317;0.36] 0.361 [0.348;0.374] 0.382 [0.364;0.402] 5.817% -5.817%
Fundão GDM - All candidate SNPs (380) 0.323 [0.307;0.343] 0.367 [0.354;0.378] 0.412 [0.39;0.433] 12.262% -11.989%
Fundão GF - All candidate SNPs (380) 0.325 [0.306;0.343] 0.366 [0.354;0.379] 0.409 [0.389;0.43] 11.749% -11.202%
Fundão LFMM - All candidate SNPs (380) 0.301 [0.281;0.321] 0.363 [0.35;0.375] 0.431 [0.409;0.454] 18.733% -17.08%
Fundão pRDA - All candidate SNPs (380) 0.312 [0.294;0.33] 0.367 [0.355;0.379] 0.426 [0.404;0.448] 16.076% -14.986%
Fundão RDA - All candidate SNPs (380) 0.3 [0.28;0.32] 0.363 [0.352;0.376] 0.431 [0.41;0.454] 18.733% -17.355%
Fundão GDM - Common candidate SNPs (69) 0.323 [0.304;0.34] 0.367 [0.355;0.379] 0.414 [0.394;0.435] 12.807% -11.989%
Fundão GF - Common candidate SNPs (69) 0.326 [0.309;0.344] 0.367 [0.354;0.379] 0.41 [0.387;0.43] 11.717% -11.172%
Fundão LFMM - Common candidate SNPs (69) 0.317 [0.298;0.338] 0.363 [0.351;0.376] 0.412 [0.391;0.434] 13.499% -12.672%
Fundão pRDA - Common candidate SNPs (69) 0.353 [0.337;0.37] 0.366 [0.354;0.378] 0.379 [0.36;0.399] 3.552% -3.552%
Fundão RDA - Common candidate SNPs (69) 0.308 [0.29;0.325] 0.367 [0.354;0.379] 0.429 [0.409;0.453] 16.894% -16.076%
Fundão GDM - Candidate SNPs considering pop. struct. correction (221) 0.314 [0.294;0.333] 0.364 [0.351;0.376] 0.418 [0.396;0.439] 14.835% -13.736%
Fundão GF - Candidate SNPs considering pop. struct. correction (221) 0.317 [0.299;0.336] 0.365 [0.352;0.377] 0.415 [0.393;0.435] 13.699% -13.151%
Fundão LFMM - Candidate SNPs considering pop. struct. correction (221) 0.304 [0.283;0.324] 0.363 [0.352;0.376] 0.427 [0.404;0.45] 17.631% -16.253%
Fundão pRDA - Candidate SNPs considering pop. struct. correction (221) 0.318 [0.301;0.337] 0.367 [0.356;0.38] 0.42 [0.398;0.443] 14.441% -13.351%
Fundão RDA - Candidate SNPs considering pop. struct. correction (221) 0.29 [0.271;0.31] 0.364 [0.352;0.376] 0.444 [0.422;0.467] 21.978% -20.33%
Fundão GDM - Control SNPs unmatching allele frequencies (380) 0.352 [0.331;0.373] 0.362 [0.35;0.375] 0.373 [0.355;0.394] 3.039% -2.762%
Fundão GF - Control SNPs unmatching allele frequencies (380) 0.319 [0.298;0.341] 0.359 [0.347;0.372] 0.402 [0.382;0.423] 11.978% -11.142%
Fundão LFMM - Control SNPs unmatching allele frequencies (380) 0.304 [0.283;0.323] 0.358 [0.347;0.37] 0.416 [0.396;0.435] 16.201% -15.084%
Fundão pRDA - Control SNPs unmatching allele frequencies (380) 0.308 [0.289;0.328] 0.367 [0.355;0.38] 0.431 [0.409;0.455] 17.439% -16.076%
Fundão RDA - Control SNPs unmatching allele frequencies (380) 0.313 [0.292;0.335] 0.357 [0.345;0.37] 0.404 [0.384;0.423] 13.165% -12.325%
Fundão GDM - Control SNPs matching allele frequencies (380) 0.35 [0.331;0.371] 0.364 [0.351;0.376] 0.377 [0.357;0.398] 3.571% -3.846%
Fundão GF - Control SNPs matching allele frequencies (380) 0.32 [0.299;0.34] 0.361 [0.348;0.373] 0.405 [0.384;0.425] 12.188% -11.357%
Fundão LFMM - Control SNPs matching allele frequencies (380) 0.3 [0.281;0.321] 0.358 [0.344;0.369] 0.42 [0.399;0.441] 17.318% -16.201%
Fundão pRDA - Control SNPs matching allele frequencies (380) 0.284 [0.264;0.303] 0.363 [0.351;0.375] 0.451 [0.428;0.473] 24.242% -21.763%
Fundão RDA - Control SNPs matching allele frequencies (380) 0.306 [0.283;0.326] 0.357 [0.345;0.37] 0.413 [0.392;0.434] 15.686% -14.286%
Fundão GDM - All SNPs (9817) 0.341 [0.322;0.362] 0.363 [0.351;0.375] 0.386 [0.366;0.407] 6.336% -6.061%
Fundão GF - All SNPs (9817) 0.316 [0.295;0.337] 0.36 [0.348;0.373] 0.407 [0.386;0.428] 13.056% -12.222%
Fundão LFMM - All SNPs (9817) 0.296 [0.276;0.318] 0.359 [0.346;0.371] 0.426 [0.406;0.448] 18.663% -17.549%
Fundão pRDA - All SNPs (9817) 0.298 [0.279;0.317] 0.366 [0.355;0.379] 0.441 [0.417;0.463] 20.492% -18.579%
Fundão RDA - All SNPs (9817) 0.303 [0.282;0.324] 0.357 [0.344;0.368] 0.414 [0.392;0.432] 15.966% -15.126%
Bordeaux CTD - Mean annual temperature (bio1, °C) 0.04 [0.032;0.047] 0.034 [0.029;0.04] 0.03 [0.022;0.037] -11.765% 17.647%
Bordeaux CTD - Annual precipitation (bio12, mm) 0.039 [0.031;0.048] 0.035 [0.03;0.04] 0.032 [0.025;0.04] -8.571% 11.429%
Bordeaux CTD - Precipitation seasonality (bio15, index) 0.03 [0.024;0.037] 0.036 [0.031;0.041] 0.043 [0.033;0.052] 19.444% -16.667%
Bordeaux CTD - Isothermality (bio3, index) 0.04 [0.03;0.049] 0.035 [0.03;0.04] 0.032 [0.024;0.039] -8.571% 14.286%
Bordeaux CTD - Temperature seasonality (bio4, °C) 0.027 [0.021;0.033] 0.035 [0.03;0.04] 0.046 [0.037;0.054] 31.429% -22.857%
Bordeaux CTD - Summer heat moisture index (SHM, °C/mm) 0.037 [0.029;0.045] 0.035 [0.03;0.04] 0.034 [0.026;0.043] -2.857% 5.714%
Bordeaux CTD - All climatic variables (Euclidean distance) 0.035 [0.028;0.042] 0.036 [0.03;0.041] 0.036 [0.027;0.046] 0% -2.778%
Bordeaux GDM - All candidate SNPs (380) 0.028 [0.022;0.033] 0.036 [0.03;0.041] 0.046 [0.038;0.055] 27.778% -22.222%
Bordeaux GF - All candidate SNPs (380) 0.028 [0.022;0.034] 0.036 [0.03;0.041] 0.045 [0.036;0.054] 25% -22.222%
Bordeaux LFMM - All candidate SNPs (380) 0.029 [0.023;0.035] 0.036 [0.031;0.042] 0.046 [0.036;0.056] 27.778% -19.444%
Bordeaux pRDA - All candidate SNPs (380) 0.033 [0.025;0.039] 0.036 [0.03;0.041] 0.039 [0.03;0.05] 8.333% -8.333%
Bordeaux RDA - All candidate SNPs (380) 0.03 [0.023;0.036] 0.037 [0.031;0.042] 0.046 [0.035;0.057] 24.324% -18.919%
Bordeaux GDM - Common candidate SNPs (69) 0.028 [0.022;0.034] 0.035 [0.03;0.041] 0.045 [0.037;0.054] 28.571% -20%
Bordeaux GF - Common candidate SNPs (69) 0.028 [0.022;0.034] 0.035 [0.03;0.04] 0.045 [0.037;0.054] 28.571% -20%
Bordeaux LFMM - Common candidate SNPs (69) 0.029 [0.023;0.035] 0.036 [0.031;0.041] 0.045 [0.037;0.054] 25% -19.444%
Bordeaux pRDA - Common candidate SNPs (69) 0.03 [0.024;0.037] 0.036 [0.031;0.042] 0.043 [0.033;0.053] 19.444% -16.667%
Bordeaux RDA - Common candidate SNPs (69) 0.036 [0.028;0.044] 0.035 [0.029;0.04] 0.035 [0.024;0.044] 0% 2.857%
Bordeaux GDM - Candidate SNPs considering pop. struct. correction (221) 0.029 [0.023;0.035] 0.036 [0.031;0.042] 0.046 [0.036;0.054] 27.778% -19.444%
Bordeaux GF - Candidate SNPs considering pop. struct. correction (221) 0.03 [0.023;0.036] 0.036 [0.03;0.041] 0.043 [0.035;0.053] 19.444% -16.667%
Bordeaux LFMM - Candidate SNPs considering pop. struct. correction (221) 0.03 [0.023;0.036] 0.036 [0.031;0.042] 0.044 [0.034;0.055] 22.222% -16.667%
Bordeaux pRDA - Candidate SNPs considering pop. struct. correction (221) 0.032 [0.026;0.039] 0.036 [0.031;0.042] 0.041 [0.03;0.051] 13.889% -11.111%
Bordeaux RDA - Candidate SNPs considering pop. struct. correction (221) 0.034 [0.026;0.042] 0.036 [0.03;0.041] 0.038 [0.026;0.049] 5.556% -5.556%
Bordeaux GDM - Control SNPs unmatching allele frequencies (380) 0.028 [0.022;0.035] 0.037 [0.031;0.042] 0.048 [0.036;0.059] 29.73% -24.324%
Bordeaux GF - Control SNPs unmatching allele frequencies (380) 0.028 [0.022;0.035] 0.036 [0.032;0.042] 0.047 [0.036;0.057] 30.556% -22.222%
Bordeaux LFMM - Control SNPs unmatching allele frequencies (380) 0.031 [0.024;0.037] 0.037 [0.031;0.042] 0.044 [0.033;0.055] 18.919% -16.216%
Bordeaux pRDA - Control SNPs unmatching allele frequencies (380) 0.035 [0.027;0.042] 0.036 [0.03;0.041] 0.037 [0.028;0.046] 2.778% -2.778%
Bordeaux RDA - Control SNPs unmatching allele frequencies (380) 0.03 [0.023;0.036] 0.036 [0.031;0.042] 0.045 [0.034;0.055] 25% -16.667%
Bordeaux GDM - Control SNPs matching allele frequencies (380) 0.028 [0.023;0.034] 0.036 [0.031;0.041] 0.047 [0.038;0.056] 30.556% -22.222%
Bordeaux GF - Control SNPs matching allele frequencies (380) 0.029 [0.023;0.035] 0.036 [0.031;0.042] 0.046 [0.035;0.055] 27.778% -19.444%
Bordeaux LFMM - Control SNPs matching allele frequencies (380) 0.031 [0.025;0.039] 0.036 [0.03;0.041] 0.042 [0.031;0.053] 16.667% -13.889%
Bordeaux pRDA - Control SNPs matching allele frequencies (380) 0.04 [0.031;0.048] 0.035 [0.029;0.04] 0.03 [0.022;0.039] -14.286% 14.286%
Bordeaux RDA - Control SNPs matching allele frequencies (380) 0.03 [0.023;0.036] 0.036 [0.031;0.042] 0.045 [0.034;0.056] 25% -16.667%
Bordeaux GDM - All SNPs (9817) 0.028 [0.022;0.034] 0.037 [0.031;0.042] 0.048 [0.038;0.058] 29.73% -24.324%
Bordeaux GF - All SNPs (9817) 0.029 [0.022;0.034] 0.037 [0.031;0.041] 0.047 [0.036;0.056] 27.027% -21.622%
Bordeaux LFMM - All SNPs (9817) 0.03 [0.023;0.037] 0.036 [0.031;0.042] 0.044 [0.033;0.054] 22.222% -16.667%
Bordeaux pRDA - All SNPs (9817) 0.036 [0.028;0.044] 0.035 [0.03;0.04] 0.035 [0.025;0.044] 0% 2.857%
Bordeaux RDA - All SNPs (9817) 0.03 [0.023;0.037] 0.036 [0.031;0.042] 0.045 [0.034;0.056] 25% -16.667%

3.5 Experimental design

We export in a table with the number and proportion of dead trees for each population in each common garden.

Code
ExpDesignTab <- lapply(unique(survival_data$site),function(site_i){
  
survival_data %>% 
    filter(site==site_i) %>% 
    dplyr::select(pop,survival) %>% 
    drop_na() %>% 
    group_by(pop) %>% 
    dplyr::summarise(nb_dead=n()-sum(survival),nb_tot=n()) %>% 
    mutate(prop_dead=nb_dead*100/nb_tot)
  
}) %>% 
  setNames(unique(survival_data$site)) %>% 
  bind_rows(.id="cg") %>% 
  pivot_wider(names_from=cg,values_from = c(nb_dead, nb_tot,prop_dead),names_sep="_") %>% 
  dplyr::select(pop,contains("asturias"),contains("bordeaux"),contains("caceres"),contains("madrid"),contains("portugal"))

# To generate a latex table
# =========================
# print(xtable(ExpDesignTab, type = "latex",digits=2), 
#       file = paste0(here("tables/ExpDesignMortalityModelsPerPopCG.tex")), 
#       include.rownames=FALSE)

# Export the table for the Supplementary Information
# =================================================
ExpDesignTab %>% 
  dplyr::select(-contains("prop_dead")) %>% 
  saveRDS(here("tables/ExpDesignMortalityModelsPerPopCG.rds"))


# Show table
# ==========
ExpDesignTab %>% kable_mydf()

# Information used in the manuscript
ExpDesignTab[,-1] %>% 
  dplyr::summarise_all(mean) %>% 
  kable_mydf

3.6 Correlation coefficients

Stan code of the mortality model without the genomic offset predictions or the CTD (i.e. only with height as predictor).

Code
stancode = stan_model(here("scripts/StanModels/ValidationCommonGarden_BinomialMortalityModel_WithoutGO.stan"))
print(stancode)
Code
lapply(unique(survival_data$site),function(site_i){
  
# Subset the data for the site i and add the initial heights
sub_data <- survival_data %>% 
  filter(site == site_i) %>% 
  group_by(pop) %>% 
  dplyr::summarise(nb_dead=n()-sum(survival),nb_tot=n()) %>% # transform survival data into mortality data
  inner_join(pop_heights %>% dplyr::select(any_of(c("height", "pop"))), by="pop") %>% 
  arrange(pop)
      
# Data in a list for Stan 
stanlist <- list(N = nrow(sub_data),
                 nb_dead = sub_data$nb_dead,
                 nb_tot = sub_data$nb_tot,
                 H = (sub_data$height-mean(sub_data$height)) / sd(sub_data$height))
    
# Running the model
mod <- sampling(stancode, data = stanlist, iter = 2000, chains = 4, cores = 4, save_warmup = FALSE) 

# Save the model and the stanlist
list(mod = mod, 
     stanlist = stanlist) %>% 
  saveRDS(file=here(paste0("outputs/ValidationCommonGarden/MortalityModels/stan_models/",
                                                           site_i,"_ModelWithoutGO.rds")))
    
})
Code
corrtab <- lapply(unique(survival_data$site),function(site_i){
  
  list_mod <- readRDS(file=here(paste0("outputs/ValidationCommonGarden/MortalityModels/stan_models/",site_i,"_ModelWithoutGO.rds")))

  # Extract posterior samples
  posterior_samples <- rstan::extract(list_mod[[1]])

  # Extract parameters (posterior samples for beta_0 and beta_H)
  beta_0_samples <- posterior_samples$beta_0
  beta_H_samples <- posterior_samples$beta_H

  # Mean tree height for each population (H)
  H <- list_mod[[2]]$H

  # Generate posterior predictions for each population
  N <- list_mod[[2]]$N # number of populations
  num_samples <- length(beta_0_samples)  # Number of posterior samples

  # Initialize matrix to store predicted probabilities
  mortality_prob <- matrix(NA, nrow=num_samples, ncol=N)

  # Compute mortality probability for each population for each posterior sample
  for (i in 1:num_samples) {
    for (j in 1:N) {
      # Logistic transformation to get probability of mortality
      mortality_prob[i, j] <- 1 / (1 + exp(-(beta_0_samples[i] + beta_H_samples[i] * H[j])))
    }
  }

  # Summary of posterior predictive mortality probabilities for each population
  mortality_prob_summary <- apply(mortality_prob, 2, function(x) c(mean=mean(x), sd=sd(x), quantile(x, probs=c(0.025, 0.975))))

  # Convert to a readable data frame (optional)
  mortality_prob_df <- tibble(
    pop = unique(df$pop),
    mean_prob = mortality_prob_summary["mean", ],
    sd_prob = mortality_prob_summary["sd", ],
    lower_95 = mortality_prob_summary["2.5%", ],
    upper_95 = mortality_prob_summary["97.5%", ],
    prob_obs = list_mod[[2]]$nb_dead / list_mod[[2]]$nb_tot, # Observed proportion of dead trees
    diff_p = prob_obs - mean_prob, # diff bwt obs prob. and estimated prob.
    cg = site_i
  ) #%>% 
    #arrange(diff_p) %>% 
    #mutate(diff_p_rank = 1:nrow(.))
  
  
  lapply(unique(df$method_input_code), function(method_input_code_i){
  
  sub_data <- df %>% 
    filter(method_input_code == method_input_code_i & cg == site_i) %>% 
    #arrange(varX) %>% 
    #mutate(varX_rank = 1:nrow(.)) %>% 
    inner_join(mortality_prob_df, by=c("pop","cg"))
  
  tibble(pearson_correlation = cor.test(sub_data$varX, sub_data$diff_p, method="pearson")$estimate[[1]],
         pearson_pvalue = cor.test(sub_data$varX, sub_data$diff_p, method="pearson")$p.value,
         spearman_correlation = cor.test(sub_data$varX, sub_data$diff_p, method="spearman")$estimate[[1]],
         spearman_pvalue = cor.test(sub_data$varX, sub_data$diff_p, method="spearman")$p.value,
         method_input_code = method_input_code_i,
         method_input_name = unique(sub_data$method_input_name),
         input_name = unique(sub_data$input_name),
         input_code = unique(sub_data$input_code),
         method = unique(sub_data$method),
         cg = site_i)
}) %>% bind_rows()
  

}) %>% bind_rows()

saveRDS(corrtab, here("outputs/ValidationCommonGarden/MortalityModels/corrtab.rds"))

We plot the correlations.

Code
corrtab <- readRDS(here("outputs/ValidationCommonGarden/MortalityModels/corrtab.rds"))

correlation_types <- c("pearson_correlation", "spearman_correlation")

lapply(correlation_types, function(coeff){
  
  p <- corrtab %>%
    pivot_longer(values_to = "correlation_estimate", names_to = "correlation_type", cols = all_of(correlation_types)) %>% 
    mutate(correlation_pvalue = ifelse(correlation_type == "pearson_correlation", pearson_pvalue, spearman_pvalue)) %>% 
    dplyr::select(-pearson_pvalue,-spearman_pvalue) %>% 
    mutate(pvalue_signi = ifelse(correlation_pvalue < 0.05, "p-value < 0.05", "p-value ≥ 0.05")) %>% 
    filter(correlation_type == coeff) %>% 
    mutate(cg_name = case_when(cg == "asturias" ~ "Asturias, Spain (37 months)",
                               cg == "bordeaux" ~ "Bordeaux, France (85 months)",
                               cg == "caceres" ~ "Cáceres, Spain (8 months)",
                               cg == "madrid" ~ "Madrid, Spain (13 months)",
                               cg == "portugal" ~ "Fundão, Portugal (27 months)")) %>% 
    mutate(input_name = factor(input_name, levels = c(unique(corrtab$input_name)[[13]],
                                                      unique(corrtab$input_name)[[11]],
                                                      unique(corrtab$input_name)[[12]],
                                                      unique(corrtab$input_name)[[10]],
                                                      unique(corrtab$input_name)[[8]],
                                                      unique(corrtab$input_name)[[9]],
                                                      unique(corrtab$input_name)[[1]],
                                                      unique(corrtab$input_name)[[2]],
                                                      unique(corrtab$input_name)[[3]],
                                                      unique(corrtab$input_name)[[4]],
                                                      unique(corrtab$input_name)[[5]],
                                                      unique(corrtab$input_name)[[6]],
                                                      unique(corrtab$input_name)[[7]])),
           method = factor(method, levels = c(unique(corrtab$method)[[2]],
                                              unique(corrtab$method)[[3]],
                                              unique(corrtab$method)[[4]],
                                              unique(corrtab$method)[[5]],
                                              unique(corrtab$method)[[6]],
                                              unique(corrtab$method)[[1]])),
           cg_name = factor(cg_name, levels = c("Madrid, Spain (13 months)",
                                                "Cáceres, Spain (8 months)",
                                                "Asturias, Spain (37 months)",
                                                "Bordeaux, France (85 months)",
                                                "Fundão, Portugal (27 months)")))

# Plots with CTD and SNP sets
# ===========================
p_allvar <- p %>% ggplot(aes(x = method,
                             y = correlation_estimate,
                             color = input_name,
                             shape = pvalue_signi)) +
  geom_rect(inherit.aes=FALSE, aes(xmin=-Inf, xmax=Inf, ymin=0, ymax=Inf), color="transparent", fill="green", alpha=0.005) + 
  geom_hline(yintercept = 0,color="gray") +
  geom_point(position = position_dodge(width = .5), size=5) +
  facet_wrap(~cg_name, ncol=2) + 
  {if(coeff=="pearson_correlation") ylab("Pearson correlation estimates")} + 
  {if(coeff=="spearman_correlation") ylab("Spearman correlation estimates (rank correlations)") } + 
  xlab("") +
  ylim(c(-0.44,0.6)) +
  scale_color_manual(values=colors_coeff)+
  scale_shape_manual(values = c(16,18)) + #, name="p-value of the correlations"
  theme_bw() +
  labs(color="",shape="") +
  theme(axis.text.x =  element_text(size=13),
        axis.text.y = element_text(size=13),
        axis.title.y = element_text(size=16),
        legend.title = element_blank(), 
        strip.text.x = element_text(size = 16),
        strip.background = element_blank(),
        legend.position = c(0.77,0.15),
        legend.text=element_text(size=13),
        panel.grid.minor.x=element_blank(),
        panel.grid.major.x=element_blank()) +
  guides(color = guide_legend(ncol=1),
         shape = guide_legend(override.aes = list(size =2 )))
p_allvar
# save in pdf and png
ggsave(p_allvar,
       file=here(paste0("figs/ValidationCommonGarden/MortalityModels/",coeff,"_SNPsandCTD.pdf")),
       device="pdf",
       height=13.6,
       width=15)

ggsave(p_allvar,
       file=here(paste0("figs/ValidationCommonGarden/MortalityModels/",coeff,"_SNPsandCTD.png")),
       height=13.6,
       width=15)



# Adding some images to represent the climatic differences among common gardens
annotation_custom2 <- function (grob, xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = Inf, data){ layer(data = data, stat = StatIdentity, position = PositionIdentity, 
        geom = ggplot2:::GeomCustomAnn,
        inherit.aes = TRUE, params = list(grob = grob, 
                                          xmin = xmin, xmax = xmax, 
                                          ymin = ymin, ymax = ymax))}

df_png <- p %>% 
  dplyr::select(cg_name) %>%
  distinct %>% 
  dplyr::mutate(image = case_when(cg_name == "Asturias, Spain (37 months)" ~ "reports/cloud.png",
                                  cg_name == "Bordeaux, France (85 months)" ~ "reports/cloud.png",
                                  cg_name == "Cáceres, Spain (8 months)" ~ "reports/sun.png",
                                  cg_name == "Madrid, Spain (13 months)" ~ "reports/sun.png",
                                  cg_name == "Fundão, Portugal (27 months)" ~ 'reports/cloud_and_sun.png'))

list_pngs <- lapply(unique(p$cg_name), function(cg_name_i){
  
sub <-  p %>% 
  filter(cg_name==cg_name_i) %>%
  slice(1)

png_cg = annotation_custom2(rasterGrob(readPNG(here(df_png[df_png$cg_name==cg_name_i,"image"])),interpolate=TRUE), 
                            ymin = -0.5,
                            ymax= -0.3,
                            xmin = 0.5,
                            xmax = 1.1, 
                            data=sub)
})

p_allvar_images <- p_allvar + list_pngs[[1]]+ list_pngs[[2]]+ list_pngs[[3]]+ list_pngs[[4]]+ list_pngs[[5]]

p_allvar_images

# save in pdf and png
ggsave(p_allvar_images,
       file=here(paste0("figs/ValidationCommonGarden/MortalityModels/",coeff,"_SNPsandCTD_WithImages.pdf")),
       device="pdf",
       height=13.6,
       width=15)

})
Warning in geom_rect(inherit.aes = FALSE, aes(xmin = -Inf, xmax = Inf, ymin = 0, : All aesthetics have length 1, but the data has 185 rows.
ℹ Please consider using `annotate()` or provide this layer with data containing a single row.
Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, : for 'p-value ≥ 0.05' in 'mbcsToSbcs': >= substituted for ≥ (U+2265)
Warning in geom_rect(inherit.aes = FALSE, aes(xmin = -Inf, xmax = Inf, ymin = 0, : All aesthetics have length 1, but the data has 185 rows.
ℹ Please consider using `annotate()` or provide this layer with data containing a single row.
All aesthetics have length 1, but the data has 185 rows.
ℹ Please consider using `annotate()` or provide this layer with data containing a single row.
Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, : for 'p-value ≥ 0.05' in 'mbcsToSbcs': >= substituted for ≥ (U+2265)
Warning in geom_rect(inherit.aes = FALSE, aes(xmin = -Inf, xmax = Inf, ymin = 0, : All aesthetics have length 1, but the data has 185 rows.
ℹ Please consider using `annotate()` or provide this layer with data containing a single row.
Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, : for 'p-value ≥ 0.05' in 'mbcsToSbcs': >= substituted for ≥ (U+2265)
Warning in geom_rect(inherit.aes = FALSE, aes(xmin = -Inf, xmax = Inf, ymin = 0, : All aesthetics have length 1, but the data has 185 rows.
ℹ Please consider using `annotate()` or provide this layer with data containing a single row.
All aesthetics have length 1, but the data has 185 rows.
ℹ Please consider using `annotate()` or provide this layer with data containing a single row.
Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, : for 'p-value ≥ 0.05' in 'mbcsToSbcs': >= substituted for ≥ (U+2265)
[[1]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/pearson_correlation_SNPsandCTD_WithImages.pdf"

[[2]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/MortalityModels/spearman_correlation_SNPsandCTD_WithImages.pdf"

4 Height models

In this section, we estimate the association between genomic offset (GO) predictions or climate transfer distances (CTD) and tree height, independently in the five common gardens. We compare four different mathematical models.

We first subset height measurements from the dataset (below the first five rows of the dataset are shown).

Code
height_measurements <- c("AST_htmar14","BDX_htnov18","CAC_htdec11","MAD_htdec11","POR_htmay13")

height_data <- pheno_data %>% 
  dplyr::rename(cg = site) %>% 
  dplyr::select(cg,block,pop,clon,tree,any_of(height_measurements)) %>% 
  pivot_longer(cols=any_of(height_measurements), names_to=NULL,values_to="height",values_drop_na = TRUE)

height_data[1:5,] %>% kable_mydf
cg block pop clon tree height
asturias 2 CAR CAR12 CAR12_2 840
asturias 4 STJ STJ14 STJ14_4 1200
asturias 5 CUE CUE6 CUE6_5 1160
asturias 6 MAD MAD1 MAD1_6 1340
asturias 4 ORI ORI12 ORI12_4 1320

We write a function to plot the model coefficients :

Code
coeff_y_axis <- list(beta_X1 = "Regression coefficients $\\beta_{X_1}$ in height models",
                     beta_X2 = "Regression coefficients $\\beta_{X_2}$ in height models",
                     beta_H = "Regression coefficients $\\beta_H$ in height models",
                     R_squared = "Proportion of variance explained ($R^2$) in height models",
                     classic_R2 = "Classic $R^2$ in height models",
                     bayes_R2_res = "Residual-based Bayes $R^2$ in height models",
                     bayes_R2_mod = "Model-based Bayes $R^2$ in height models")


generate_interval_plots <- function(model_i){
  
  coefftab <- readRDS(file=here(paste0("outputs/ValidationCommonGarden/HeightModels/coefftab_",model_i,".rds")))
  params <-  coefftab$term %>% unique() %>% str_subset("^beta|squared|R2")
  
  
  p <- lapply(params, function(coeff){
    
    if(coeff %in% c("beta_X1", "beta_X2","beta_H","R_squared")){
      p_beta <- coefftab %>% 
        filter(term %in% c("beta_X1", "beta_X2")) %>% 
        pivot_wider(names_from = "term", values_from = c("mean","conf_low","conf_high","std_deviation")) %>% 
        mutate(expected = ifelse(conf_high_beta_X1<0 & conf_low_beta_X2<0 & conf_high_beta_X2>0 , "Expected associations", "Non-expected associations")) %>% 
        dplyr::select(-contains("mean"),-contains("conf"), -contains("deviation"))
      }
    
    p <- coefftab %>% filter(term==coeff) 
    
    if(coeff %in% c("beta_X1", "beta_X2","beta_H","R_squared")){p <- p %>% left_join(p_beta)}
    
    p <- p %>% mutate(cg_name = case_when(cg == "asturias" ~ "Asturias, Spain (37 months)",
                                          cg == "bordeaux" ~ "Bordeaux, France (85 months)",
                                          cg == "caceres" ~ "Cáceres, Spain (8 months)",
                                          cg == "madrid" ~ "Madrid, Spain (13 months)",
                                          cg == "portugal" ~ "Fundão, Portugal (27 months)"),
                      input_name = factor(input_name, levels = c(unique(coefftab$input_name)[[13]],
                                                                 unique(coefftab$input_name)[[11]],
                                                                 unique(coefftab$input_name)[[12]],
                                                                 unique(coefftab$input_name)[[10]],
                                                                 unique(coefftab$input_name)[[8]],
                                                                 unique(coefftab$input_name)[[9]],
                                                                 unique(coefftab$input_name)[[1]],
                                                                 unique(coefftab$input_name)[[2]],
                                                                 unique(coefftab$input_name)[[3]],
                                                                 unique(coefftab$input_name)[[4]],
                                                                 unique(coefftab$input_name)[[5]],
                                                                 unique(coefftab$input_name)[[6]],
                                                                 unique(coefftab$input_name)[[7]])),
                      method = factor(method, levels = c(unique(coefftab$method)[[2]],
                                                         unique(coefftab$method)[[3]],
                                                         unique(coefftab$method)[[4]],
                                                         unique(coefftab$method)[[5]],
                                                         unique(coefftab$method)[[6]],
                                                         unique(coefftab$method)[[1]])),
                      cg_name = factor(cg_name, levels = c("Madrid, Spain (13 months)",
                                                           "Cáceres, Spain (8 months)",
                                                           "Asturias, Spain (37 months)",
                                                           "Bordeaux, France (85 months)",
                                                           "Fundão, Portugal (27 months)")))

# Plots with CTD and SNP sets
# ===========================
    if(coeff %in% c("beta_X1", "beta_X2","beta_H", "R_squared")){
      p_allvar <- p %>% 
        ggplot(aes(x = method,
                   y = mean,
                   ymin = conf_low, 
                   ymax = conf_high,
                   color = input_name,
                   shape = expected)) +
        {if(coeff=="beta_X1")  geom_rect(inherit.aes=FALSE, aes(xmin=-Inf, xmax=Inf, ymin=-Inf, ymax=0), color="transparent", fill="green", alpha=0.005)} +
        {if(coeff %in% c("classic_R2", "bayes_R2_mod", "bayes_R2_res") & model_i=="M2")  geom_rect(data=nullR2,inherit.aes=FALSE, aes(xmin=-Inf, xmax=Inf,
                                                                                                                                      ymin=conf_low, ymax=conf_high),
                                                                                                   color="transparent", fill="brown", alpha=0.2)} + 
        {if(coeff %in% c("classic_R2", "bayes_R2_mod", "bayes_R2_res") & model_i=="M2")  geom_hline(data=nullR2,aes(yintercept=mean), color="brown")} +
        geom_hline(yintercept = 0,color="gray") +
        geom_pointinterval(position = position_dodge(width = .4),point_size=3.5,size=3) + # 
        facet_wrap(~cg_name, ncol=2) + #, scales="free_x", space = "free" 
        ylab(TeX(coeff_y_axis[[coeff]])) + xlab("") +
        scale_color_manual(values=colors_coeff) +
        scale_shape_manual(
          values = c(16, 3),
          labels = c(expression("Expected associations, i.e., " ~ CI[95~"%"](beta[X1]) < 0 ~ "& " ~ 0 %in% CI[95~"%"](beta[X2])),
                     "Non-expected associations")) +
        theme_bw() +
        labs(color="",shape="") +
        theme(axis.text.x =  element_text(size=13),
              axis.text.y = element_text(size=13),
              axis.title.y = element_text(size=16),
              axis.title.x = element_text(size=1),
              legend.title=element_blank(), 
              strip.text.x = element_text(size = 16),
              strip.background = element_blank(),
              legend.position = c(0.77,0.15),
              legend.box.spacing = unit(0, "pt"),
              legend.text=element_text(size=11),
              panel.grid.minor.x=element_blank(),
              panel.grid.major.x=element_blank()) +
        guides(color=guide_legend(ncol=1),
               shape = guide_legend(override.aes = list(size =2)))
      
       # save in pdf
      ggsave(p_allvar,
             file=here(paste0("figs/ValidationCommonGarden/HeightModels/",coeff,"_",model_i,"_SNPsandCTD.pdf")),
             device="pdf",
             height=13.6,
             width=15)
      
      
      
      # Figure in the main manuscript
      ###############################
      
      # We want to add some images to represent the climatic differences among common gardens
      
      
      if(coeff=="beta_X1"& model_i=="M2"){
        
        annotation_custom2 <- function (grob, xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = Inf, data){ layer(data = data, 
                                                                                                             stat = StatIdentity, 
                                                                                                             position = PositionIdentity, 
                                                                                                             geom = ggplot2:::GeomCustomAnn,
                                                                                                             inherit.aes = TRUE, 
                                                                                                             params = list(grob = grob, 
                                                                                                                           xmin = xmin, 
                                                                                                                           xmax = xmax, 
                                                                                                                           ymin = ymin, 
                                                                                                                           ymax = ymax))}
        
        df_png <- p %>% 
          dplyr::select(cg_name) %>%
          distinct %>% 
          dplyr::mutate(image = case_when(cg_name == "Asturias, Spain (37 months)" ~ "reports/cloud.png",
                                          cg_name == "Bordeaux, France (85 months)" ~ "reports/cloud.png",
                                          cg_name == "Cáceres, Spain (8 months)" ~ "reports/sun.png",
                                          cg_name == "Madrid, Spain (13 months)" ~ "reports/sun.png",
                                          cg_name == "Fundão, Portugal (27 months)" ~ 'reports/cloud_and_sun.png'))
        
        list_pngs <- lapply(unique(p$cg_name), function(cg_name_i){
          sub <-  p %>% 
            filter(cg_name==cg_name_i) %>%
            slice(1)
          
          png_cg = annotation_custom2(rasterGrob(readPNG(here(df_png[df_png$cg_name==cg_name_i,"image"])),interpolate=TRUE), 
                                      ymin = -0.44,
                                      ymax= -0.3,
                                      xmin = 0.37,
                                      xmax = 1.4, 
                                      data=sub)
          })
        
        
        p_allvar_images <- p_allvar + list_pngs[[1]]+ list_pngs[[2]]+ list_pngs[[3]]+ list_pngs[[4]]+ list_pngs[[5]]
        
        # save in pdf
        ggsave(p_allvar_images,
               file=here(paste0("figs/ValidationCommonGarden/HeightModels/",coeff,"_",model_i,"_SNPsandCTD_WithCloudAndSunImages.pdf")),
               device="pdf",
               height=13.6,
               width=15)
      }
      }
    
    if(coeff %in% c("classic_R2", "bayes_R2_mod", "bayes_R2_res") & model_i=="M2"){
      nullR2 <- readRDS(file = here("outputs/ValidationCommonGarden/HeightModels/coefftab_M2_WithoutPredictor.rds")) %>%
        filter(term==coeff) %>% 
        inner_join(distinct(p[,c("cg","cg_name")]),by="cg")
      
      
      p_allvar <- p %>% 
        ggplot(aes(x = method,
                   y = mean,
                   ymin = conf_low, 
                   ymax = conf_high,
                   color = input_name,
                   shape = input_name)) +
        {if(coeff=="beta_X1")  geom_rect(inherit.aes=FALSE, aes(xmin=-Inf, xmax=Inf, ymin=-Inf, ymax=0), color="transparent", fill="green", alpha=0.005)} + 
        {if(coeff %in% c("classic_R2", "bayes_R2_mod", "bayes_R2_res") & model_i=="M2")  geom_rect(data=nullR2,inherit.aes=FALSE, aes(xmin=-Inf, xmax=Inf,
                                                                                                                                      ymin=conf_low, ymax=conf_high),
                                                                                                   color="transparent", fill="brown", alpha=0.2)} + 
        {if(coeff %in% c("classic_R2", "bayes_R2_mod", "bayes_R2_res") & model_i=="M2")  geom_hline(data=nullR2,aes(yintercept=mean), color="brown")} +
        geom_hline(yintercept = 0,color="gray") +
        geom_pointinterval(position = position_dodge(width = .4),point_size=3.5,size=3) + # 
        facet_wrap(~cg_name, ncol=2) + #, scales="free_x", space = "free" 
        ylab(TeX(coeff_y_axis[[coeff]])) + xlab("") +
        scale_color_manual(values=colors_coeff) +
        scale_shape_manual(values = c(rep(17,6),rep(16,7))) +
        theme_bw() +
        labs(color="",shape="") +
        theme(axis.text.x =  element_text(size=13),
              axis.text.y = element_text(size=13),
              axis.title.y = element_text(size=16),
              axis.title.x = element_text(size=1),
              legend.title=element_blank(), 
              strip.text.x = element_text(size = 16),
              strip.background = element_blank(),
              legend.position = c(0.77,0.15),
              legend.text=element_text(size=11),
              panel.grid.minor.x=element_blank(),
              panel.grid.major.x=element_blank()) +
        guides(color=guide_legend(ncol=1),
               shape = guide_legend(override.aes = list(size =2 )))
      
          # save in pdf
      ggsave(p_allvar,
             file=here(paste0("figs/ValidationCommonGarden/HeightModels/",coeff,"_",model_i,"_SNPsandCTD.pdf")),
             device="pdf",
             height=11,
             width=11)
}
    

    
    # save in png
    # ggsave(p_allvar,
    #        file=here(paste0("figs/ValidationCommonGarden/HeightModels/",coeff,"_",model_i,"_SNPsandCTD.png")),
    #        height=11,
    #        width=14)
    
    
    # Plots only SNP sets
    # ===================
    if(coeff %in% c("beta_X1", "beta_X2", "beta_H",  "R_squared")){

            p_SNPs <- p %>% 
              filter(!method == "CTD") %>% 
              ggplot(aes(x = method,
                         y = mean,
                         ymin = conf_low, 
                         ymax = conf_high,
                         color = input_name,
                         shape = expected)) +
              {if(coeff=="beta_X1")  geom_rect(inherit.aes=FALSE, aes(xmin=-Inf, xmax=Inf, ymin=-Inf, ymax=0), color="transparent", fill="green", alpha=0.005)} +
              {if(coeff %in% c("classic_R2", "bayes_R2_mod", "bayes_R2_res") & model_i=="M2")  geom_rect(data=nullR2,inherit.aes=FALSE, aes(xmin=-Inf, 
                                                                                                                                            xmax=Inf,
                                                                                                                                            ymin=conf_low,
                                                                                                                                            ymax=conf_high),
                                                                                                         color="transparent", fill="brown", alpha=0.2)} + 

              {if(coeff %in% c("classic_R2", "bayes_R2_mod", "bayes_R2_res") & model_i=="M2")  geom_hline(data=nullR2,aes(yintercept=mean), color="brown")} + 
              geom_hline(yintercept = 0,color="gray") +
              geom_pointinterval(position = position_dodge(width = .4),point_size=3.5,size=3) +
              facet_wrap(~cg_name, ncol=2) + #, scales="free_x", space = "free" 
              ylab(TeX(coeff_y_axis[[coeff]])) + xlab("") +
              scale_color_manual(values=colors_coeff) +
              scale_shape_manual(
                values = c(16, 3),
                labels = c(expression("Expected associations, i.e., " ~ CI[95~"%"](beta[X1]) < 0 ~ "& " ~ 0 %in% CI[95~"%"](beta[X2])),
                           "Non-expected associations"))+
              theme_bw() +
              labs(color="",shape="") +
              theme(axis.text.x =  element_text(size=13),
                    axis.text.y = element_text(size=13),
                    axis.title.y = element_text(size=16),
                    axis.title.x = element_text(size=1),
                    legend.title=element_blank(), 
                    strip.text.x = element_text(size = 16),
                    strip.background = element_blank(),
                    legend.position = c(0.77,0.15),
                    legend.text=element_text(size=13),
                    panel.grid.minor.x=element_blank(),
                    panel.grid.major.x=element_blank()) +
              guides(color=guide_legend(ncol=1),
                     shape = guide_legend(override.aes = list(size =2 )))
            
                # save in pdf
            ggsave(p_SNPs,
                   file=here(paste0("figs/ValidationCommonGarden/HeightModels/",coeff,"_",model_i,"_SNPsets.pdf")),
                   device="pdf",
                   height=12,
                   width=12)
    
    }
    
    
    
    if(coeff %in% c("classic_R2", "bayes_R2_mod", "bayes_R2_res") & model_i=="M2"){
      
      p_SNPs <- p %>% 
        filter(!method == "CTD") %>% 
        ggplot(aes(x = method,
                   y = mean,
                   ymin = conf_low, 
                   ymax = conf_high,
                   color = input_name,
                   shape = input_name)) +
        {if(coeff=="beta_X1")  geom_rect(inherit.aes=FALSE, aes(xmin=-Inf, xmax=Inf, ymin=-Inf, ymax=0), color="transparent", fill="green", alpha=0.005)} + 
        {if(coeff %in% c("classic_R2", "bayes_R2_mod", "bayes_R2_res") & model_i=="M2")  geom_rect(data=nullR2,inherit.aes=FALSE, aes(xmin=-Inf, 
                                                                                                                                      xmax=Inf,
                                                                                                                                      ymin=conf_low, 
                                                                                                                                      ymax=conf_high),
                                                                                                   color="transparent", fill="brown", alpha=0.2)} + 
        {if(coeff %in% c("classic_R2", "bayes_R2_mod", "bayes_R2_res") & model_i=="M2")  geom_hline(data=nullR2,aes(yintercept=mean), color="brown")} +
        geom_hline(yintercept = 0,color="gray") +
        geom_pointinterval(position = position_dodge(width = .4),point_size=3.5,size=3) +
        facet_wrap(~cg_name, ncol=2) + #, scales="free_x", space = "free" 
        ylab(TeX(coeff_y_axis[[coeff]])) + xlab("") +
        scale_color_manual(values=colors_coeff) +
        scale_shape_manual(values = c(rep(17,6),rep(16,7))) +
        theme_bw() +
        labs(color="",shape="") +
        theme(axis.text.x =  element_text(size=13),
              axis.text.y = element_text(size=13),
              axis.title.y = element_text(size=16),
              axis.title.x = element_text(size=1),
              legend.title=element_blank(), 
              strip.text.x = element_text(size = 16),
              strip.background = element_blank(),
              legend.position = c(0.77,0.15),
              legend.text=element_text(size=13),
              panel.grid.minor.x=element_blank(),
              panel.grid.major.x=element_blank()) +
        guides(color=guide_legend(ncol=1),
               shape = guide_legend(override.aes = list(size =2 )))
      
      
      # save in pdf
      ggsave(p_SNPs,
             file=here(paste0("figs/ValidationCommonGarden/HeightModels/",coeff,"_",model_i,"_SNPsets.pdf")),
             device="pdf",
             height=11,
             width=11)
      }

    
    # save in png
    # ggsave(p_SNPs,
    #        file=here(paste0("figs/ValidationCommonGarden/HeightModels/",coeff,"_",model_i,"_SNPsets.png")),
    #        height=7,
    #        width=8)
    
    
    return(p_allvar)
    
    })
  
  return(p)
  }

We write a function to plot the quadratic relationships between tree height and the GO predictions / CTD.

Code
generate_poly_plots <- function(model_i){

coefftab <- readRDS(file=here(paste0("outputs/ValidationCommonGarden/HeightModels/coefftab_",model_i,".rds")))

x <- seq(-2,2,0.01)
poly_function <- function(a,b,x) {a * x^2 + b * x}

ggtab <- coefftab %>% 
  filter(term %in% c("beta_X1","beta_X2")) %>% 
  mutate(cg_name = case_when(cg == "asturias" ~ "Asturias, Spain (37 months)",
                             cg == "bordeaux" ~ "Bordeaux, France (85 months)",
                             cg == "caceres" ~ "Cáceres, Spain (8 months)",
                             cg == "madrid" ~ "Madrid, Spain (13 months)",
                             cg == "portugal" ~ "Fundão, Portugal (27 months)"),
         input_name = factor(input_name, levels = c(unique(coefftab$input_name)[[13]],
                                                    unique(coefftab$input_name)[[11]],
                                                    unique(coefftab$input_name)[[12]],
                                                    unique(coefftab$input_name)[[10]],
                                                    unique(coefftab$input_name)[[8]],
                                                    unique(coefftab$input_name)[[9]],
                                                    unique(coefftab$input_name)[[1]],
                                                    unique(coefftab$input_name)[[2]],
                                                    unique(coefftab$input_name)[[3]],
                                                    unique(coefftab$input_name)[[4]],
                                                    unique(coefftab$input_name)[[5]],
                                                    unique(coefftab$input_name)[[6]],
                                                    unique(coefftab$input_name)[[7]])),
         method = factor(method, levels = c(unique(coefftab$method)[[2]],
                                            unique(coefftab$method)[[3]],
                                            unique(coefftab$method)[[4]],
                                            unique(coefftab$method)[[5]],
                                            unique(coefftab$method)[[6]],
                                            unique(coefftab$method)[[1]])),
         cg_name = factor(cg_name, levels = c("Madrid, Spain (13 months)",
                                              "Cáceres, Spain (8 months)",
                                              "Asturias, Spain (37 months)",
                                              "Bordeaux, France (85 months)",
                                              "Fundão, Portugal (27 months)"))) %>% 
  pivot_wider(names_from="term", values_from = c("mean","std_deviation","conf_low","conf_high"))

ggtab <- lapply(unique(ggtab$cg_name), function(cg_name_i){
  
  lapply(unique(ggtab$method_input_name), function(method_input_name_i){
  
  subggtab <- ggtab %>% filter(method_input_name == method_input_name_i & cg_name == cg_name_i)
  
  poly_function(a=subggtab$mean_beta_X2, b=subggtab$mean_beta_X1,x=x)
  }) %>% 
    setNames(unique(ggtab$method_input_name)) %>% 
    bind_rows(.id="method_input_name") %>% 
    mutate(x_axis = x)
}) %>% 
  setNames(unique(ggtab$cg_name)) %>% 
  bind_rows(.id="cg_name") %>% 
  pivot_longer(cols=any_of(unique(ggtab$method_input_name)),names_to="method_input_name",values_to="predictions") %>% 
  inner_join(ggtab[,c("method","input_name", "method_input_name","cg","cg_name")] %>% distinct(), by=c("method_input_name","cg_name"))
  

# Below, you can use unique(ggtab$cg) or unique(ggtab$cg_legend)
ggplots <- lapply(unique(ggtab$cg_name), function(cg_name_i){
  
  lapply(unique(ggtab$method), function(method_i){
    
    ptab <- ggtab %>%
      filter(cg_name %in% cg_name_i & 
             method %in% method_i)
    
 p <- ptab %>%  ggplot(aes(x = x_axis,
             y = predictions,
             col = input_name)) +
      geom_hline(yintercept=0,color="black") +
      geom_line(linewidth=1) +
      ylab('Standardized tree height') + 
      ylim(c(min(ggtab$predictions),max(ggtab$predictions))) +
      
      {if(method_i %in% c("GDM", "GF", "RDA", "pRDA", "LFMM")) scale_color_manual(values=c("#004586FF","#FF7F0FFF","#FFD320FF","#579D1CFF","#7E0021FF","#83CAFFFF"), name = "SNP set")} +
         
      {if(method_i %in% c("GDM", "GF", "RDA", "pRDA", "LFMM")) xlab(TeX("Standardized genomic offset predictions"))} + 
      {if(method_i == "CTD") scale_color_manual(values=c("aquamarine3","#FEB5A2FF","#F02720FF","darkviolet","deeppink","cyan2","chocolate4"), name = "Climatic distance")} +
      {if(method_i == "CTD") xlab(TeX("Standardized climatic transfer distance"))} + 
      labs(title=(paste0(cg_name_i, " - ", method_i))) +
      theme_bw() +
      theme(axis.text.x = element_text(size=13),
            axis.text.y = element_text(size=13),
            axis.title = element_text(size=16),
            legend.box="horizontal",
            legend.key.width = unit(1,"cm"),
            legend.background = element_blank(),
            legend.box.background = element_blank(),
            legend.key = element_blank(),
            legend.position = c(0.68,0.83)) + 
      guides(color=guide_legend(ncol=1))

    ggsave(filename = here(paste0("figs/ValidationCommonGarden/HeightModels/ScatterPlotsPredictedRelationships/ScatterPlotsPredictedRelationships_",
                              model_i,"_",unique(ptab$cg),"_",method_i,".pdf")), device="pdf", width=7,height=7)    
  })
})
}

4.1 Model 1

4.1.1 Model equation and code

\[\begin{align*} Y_{ipb} &\sim \mathcal{N}(\mu_{pb},\sigma^{2}) \\ \mu_{pb} &= B_b + \beta_{X1}X_p + \beta_{X2}X^2_p\\ \sigma &\sim \text{Exp}(1) \\ \begin{bmatrix} B_b \\ \beta_{X1} \\ \beta_{X2} \end{bmatrix} &\sim \mathcal{N}(0,1) \\ \end{align*}\]

  • \(Y_{ipb}\) is the height of the individual \(i\) in the population \(p\) and block \(b\).
  • \(\sigma^{2}\) is the residual variance of the model.
  • \(B_b\) are the block intercepts.
  • \(X_p\) is the GO or CTD of the population \(p\), with \(\beta_{X1}\) and \(\beta_{X2}\) being its linear and quadratic coefficients, respectively. The quadratic term was included to allow for potential nonlinearity in the response, following Fitzpatrick et al. (2021).
Code
stancode = stan_model(here("scripts/StanModels/ValidationCommonGarden_HeightModel_M1.stan"))
print(stancode)
S4 class stanmodel 'anon_model' coded as follows:
data {
  int N;                                    // Number of individuals
  vector[N] Y;                              // Response variable (individual tree height)
  vector[N] X;                              // Genomic offset or climatic transfer distances
  int<lower=0> nb_bloc;                     // Number of blocks
  int<lower=0, upper=nb_bloc> bloc[N];      // Blocks
}
parameters {
  vector[nb_bloc] alpha_bloc;               // Intercepts of the blocks
  real beta_X1;                             // Linear coefficent of GO or CTD
  real beta_X2;                             // Quadratic coefficent of GO or CTD
  real<lower = 0>  sigma;                   // Residual variance of the model
}
transformed parameters {
  vector[N] mu;                             // Linear predictor
  real R_squared;                           // R^2 to evaluate the goodness of fit of the model

  mu = alpha_bloc[bloc] + beta_X1 * X + beta_X2 * square(X);
  R_squared = 1 - variance(Y - mu) / variance(Y);
}
model {

  // Likelihood
  Y ~ normal(mu, sigma);

  // Priors
  sigma ~ exponential(1);
  alpha_bloc ~ std_normal();
  beta_X1 ~ std_normal();
  beta_X2 ~ std_normal();
}
// generated quantities{
//   // log likelihood for loo
//   vector[N] log_lik;
//   vector[N] muhat;
//   for (n in 1:N) {
//     log_lik[n] = normal_lpdf(Y[n] |mu[n],sigma);
//     muhat[n] = normal_rng(mu[n], sigma);
//   }
// } 
Code
params_to_estimate <- c("beta_X1","beta_X2","R_squared","sigma","alpha_bloc")

4.1.2 Run the models

Code
lapply(unique(height_data$cg), function(site_i){
  
  lapply(unique(df$method_input_code), function(method_input_code_i){
  
    # Subset the data
    df_sub <- df %>% filter(method_input_code == method_input_code_i & cg == site_i)
    
    sub_data <- height_data %>% 
      filter(cg == site_i) %>% 
      left_join(df_sub, by = c("cg","pop"))
      
    # Data in a list for Stan
    stanlist <- list(N = nrow(sub_data),
                     Y=(sub_data$height-mean(sub_data$height))/sd(sub_data$height),
                     X=(sub_data$varX -mean(sub_data$varX))/sd(sub_data$varX),
                     nb_bloc = length(unique(sub_data$block)),
                     bloc = as.numeric(as.factor(sub_data$block)))

    # Run the models
    mod <- sampling(stancode, 
                    data = stanlist, 
                    iter = 2000, 
                    warmup = 1400, 
                    chains = 4, 
                    cores = 4, 
                    save_warmup = FALSE,
                    pars = params_to_estimate) 
    
    # Save the model and the stanlist
    list(mod = mod, 
         stanlist = stanlist) %>% 
      saveRDS(file=here(paste0("outputs/ValidationCommonGarden/HeightModels/stan_models/M1/",site_i,"_",method_input_code_i,".rds")))
    
    
  })
})
Code
coefftab <- lapply(unique(height_data$cg),function(site_i){
  
  lapply(unique(df$method_input_code), function(method_input_code_i){
    
    # Subset the data - keeping only one set of method / SNP set
    sub_data <- df %>% filter(method_input_code == method_input_code_i & cg == site_i)
  
    mod <- readRDS(file=here(paste0("outputs/ValidationCommonGarden/HeightModels/stan_models/M1/",
                                    site_i,"_",method_input_code_i,".rds")))[["mod"]]
                     
    # Save coefficients
    broom.mixed::tidyMCMC(mod,
                          droppars = NULL, 
                          robust = FALSE, # give mean and standard deviation
                          ess = F, 
                          rhat = F, 
                          conf.int = T,
                          conf.level = 0.95) %>% 
      #filter(str_detect(term, c('beta'))) %>% 
      dplyr::rename(mean=estimate,
                    std_deviation=std.error,
                    conf_low=conf.low,
                    conf_high=conf.high) %>% 
      mutate(cg = site_i,
             method_input_code = method_input_code_i,
             method_input_name = unique(sub_data$method_input_name),
             method = unique(sub_data$method),
             input_name = unique(sub_data$input_name),
             input_code = unique(sub_data$input_code),
             .before=1)     

}) %>% bind_rows()
 
}) %>% bind_rows()
  
coefftab %>% 
  saveRDS(file=here("outputs/ValidationCommonGarden/HeightModels/coefftab_M1.rds"))

4.1.3 Model coefficients

We plot the mean and 95% credible intervals of the \(\beta_{X_1}\) and \(\beta_{X_2}\) coefficients. Graph titles include the time in months corresponding to the age at which height and survival were recorded. Coefficients in the green area have the expected sign, reflecting lower height in populations with higher GO predictions.

Code
generate_interval_plots(model_i = "M1")
Joining with `by = join_by(cg, method_input_code, method_input_name, method, input_name, input_code)`
Warning in geom_rect(inherit.aes = FALSE, aes(xmin = -Inf, xmax = Inf, ymin = -Inf, : All aesthetics have length 1, but the data has 185 rows.
ℹ Please consider using `annotate()` or provide this layer with data containing a single row.
Warning in geom_rect(inherit.aes = FALSE, aes(xmin = -Inf, xmax = Inf, ymin = -Inf, : All aesthetics have length 1, but the data has 150 rows.
ℹ Please consider using `annotate()` or provide this layer with data containing a single row.
Joining with `by = join_by(cg, method_input_code, method_input_name, method, input_name, input_code)`
Joining with `by = join_by(cg, method_input_code, method_input_name, method, input_name, input_code)`
[[1]]
Warning in geom_rect(inherit.aes = FALSE, aes(xmin = -Inf, xmax = Inf, ymin = -Inf, : All aesthetics have length 1, but the data has 185 rows.
ℹ Please consider using `annotate()` or provide this layer with data containing a single row.


[[2]]


[[3]]

4.1.4 Predicted quadratic relationships

We plot the predicted quadratic relationships between tree height and GO predictions / CTD.

Code
generate_poly_plots(model_i="M1")

4.1.5 Interpretation

Higher genomic offset predictions are consistently associated with lower tree height only in Asturias. Most genomic offset predictions are also associated with lower tree height in Bordeaux.

In Fundão, Cáceres and Madrid, the association between genomic offset predictions and tree height are often in the opposite direction as expected.

Interestingly, the climatic transfer distances based in the mean annual temperature and the annual precipitation are negatively associated with tree height in all common gardens (except Cáceres for the CTD based on annual precipitation). These CTD may therefore be better predictors of tree height in common gardens than the genomic offset predictions.

Finally, I think that using tree height as a proxy for fitness to evaluate genomic offset predictions in common gardens may not be appropriate in maritime pine because this trait has a very low genetic-by-environment interaction (see previous papers). Therefore, the association between genomic offset predictions and tree height in Asturias and Bordeaux may not be due to a higher fitness of the trees that are best adapted to the climate in these common gardens, but more probably to height differences that are common across all environments (i.e. populations from Atlantic climates are generally taller). This is confirmed by the models below in which the population-specific BLUPs for height across all common gardens was included as a confounder. When included, the association between genomic offset predictions and tree height in Asturias almost disappears and either disappears or goes in the opposite direction as expected in Bordeaux. Except the genomic offset predictions based on the RDA remain negatively associated with tree height in Bordeaux and Asturias.

Therefore, I am not sure we can go very far in this validation step. We may retain that genomic offset predictions seem to show the most consistent association with tree height when generated with the RDA, and that this association is robust (though smaller) even when the fixed genetic differences across populations are included as confounder.

4.2 Model 2

4.2.1 Model equation and code

Model 2 = Model 1 + Initial height as a confounder.

\[\begin{align*} Y_{ipb} &\sim \mathcal{N}(\mu_{pb},\sigma^{2}) \\ \mu_{pb} &= B_b + \beta_{X1}X_p + \beta_{X2}X^2_p + \beta_H H_p\\ \sigma &\sim \text{Exp}(1) \\ \begin{bmatrix} B_b \\ \beta_{X1} \\ \beta_{X2} \\ \beta_H \end{bmatrix} &\sim \mathcal{N}(0,1) \\ \end{align*}\]

  • \(Y_{ipb}\) is the height of the individual \(i\) in the population \(p\) and block \(b\).
  • \(\sigma^{2}\) is the residual variance of the model.
  • \(B_b\) are the block intercepts.
  • \(X_p\) is the GO or CTD of the population \(p\), with \(\beta_{X1}\) and \(\beta_{X2}\) being its linear and quadratic coefficients, respectively. The quadratic term was included to allow for potential nonlinearity in the response, following Fitzpatrick et al. (2021).
  • \(H_p\) is a proxy of the initial height of the trees from population \(p\), i.e. when trees were planted. For that, we used the population varying intercepts calculated across all common gardens in the model 1 of Archambeau et al. (2022).
Code
stancode = stan_model(here("scripts/StanModels/ValidationCommonGarden_HeightModel_M2.stan"))
print(stancode)
S4 class stanmodel 'anon_model' coded as follows:
data {
  int N;
  vector[N] Y;                              // Response variable (individual height)
  vector[N] X;                              // Genomic offset or climatic transfer distances
  vector[N] H;                              // Initial height of the populations
  int<lower=0> nb_bloc;                     // Number of blocks
  int<lower=0, upper=nb_bloc> bloc[N];      // Blocks
}
parameters {
  vector[nb_bloc] alpha_bloc;               // Intercepts of the blocks
  real beta_X1;                             // Linear coefficent of GO or CTD
  real beta_X2;                             // Quadratic coefficent of GO or CTD
  real beta_H;                              // Coefficient of the initial height of the populations
  real<lower = 0>  sigma;                   // Residual variance of the model

}
transformed parameters {
  vector[N] mu;     // Linear predictor
  real classic_R2;  // classic R2 to evaluate the goodness of fit of the model

  mu = alpha_bloc[bloc] + beta_X1 * X + beta_X2 * square(X) + beta_H * H;
  classic_R2 = 1 - variance(Y - mu) / variance(Y); // classical R2
}
model {

  // Likelihood
  Y ~ normal(mu, sigma);

  // Priors
  sigma ~ exponential(1);
  alpha_bloc ~ std_normal();
  beta_H ~ std_normal();
  beta_X1 ~ std_normal();
  beta_X2 ~ std_normal();
}
generated quantities{
  // Bayesian R2
  real bayes_R2_res; // Residual-based based R2 (uses draws from the residual distribution)
  real bayes_R2_mod; // Model-based R2 (uses draws from the modeled residual variances)

  bayes_R2_res = variance(mu) / (variance(mu) + variance(Y-mu));
  bayes_R2_mod = variance(mu) / (variance(mu) + square(sigma));

//   // log likelihood for loo
//   vector[N] log_lik;
//   vector[N] muhat;
//   for (n in 1:N) {
//     log_lik[n] = normal_lpdf(Y[n] |mu[n],sigma);
//     muhat[n] = normal_rng(mu[n], sigma);
//   }
} 
Code
params_to_estimate <- c("beta_X1","beta_X2","beta_H","sigma","alpha_bloc",
                        "classic_R2","bayes_R2_mod","bayes_R2_res")

4.2.2 Run the models

Code
lapply(unique(height_data$cg), function(site_i){
  
  lapply(unique(df$method_input_code), function(method_input_code_i){
  
    # Subset the data
    df_sub <- df %>% filter(method_input_code == method_input_code_i & cg == site_i)
    
    sub_data <- height_data %>% 
      filter(cg == site_i) %>% 
      left_join(df_sub, by = c("cg","pop")) %>% 
      left_join(pop_heights %>% dplyr::select(any_of(c("height", "pop"))) %>% dplyr::rename(init_height=height), by="pop")
      
    # Data in a list for Stan
    stanlist <- list(N = nrow(sub_data),
                     Y=(sub_data$height-mean(sub_data$height))/sd(sub_data$height),
                     X=(sub_data$varX -mean(sub_data$varX))/sd(sub_data$varX),
                     H=(sub_data$init_height -mean(sub_data$init_height))/sd(sub_data$init_height),
                     nb_bloc = length(unique(sub_data$block)),
                     bloc = as.numeric(as.factor(sub_data$block)))

    # Run the models
    mod <- sampling(stancode, 
                    data = stanlist, 
                    iter = 2000, 
                    warmup = 1400, 
                    chains = 4, 
                    cores = 4, 
                    save_warmup = FALSE,
                    pars = params_to_estimate) 
    
    # Save the model and the stanlist
    list(mod = mod, 
         stanlist = stanlist) %>% 
      saveRDS(file=here(paste0("outputs/ValidationCommonGarden/HeightModels/stan_models/M2/",site_i,"_",method_input_code_i,".rds")))
    
    
  })
})
Code
coefftab <- lapply(unique(height_data$cg),function(site_i){
  
  lapply(unique(df$method_input_code), function(method_input_code_i){
    
    # Subset the data - keeping only one set of method / SNP set
    sub_data <- df %>% filter(method_input_code == method_input_code_i & cg == site_i)
  
    mod <- readRDS(file=here(paste0("outputs/ValidationCommonGarden/HeightModels/stan_models/M2/",
                                    site_i,"_",method_input_code_i,".rds")))[["mod"]]
                     
    # Save coefficients
    broom.mixed::tidyMCMC(mod,
                          droppars = NULL, 
                          robust = FALSE, # give mean and standard deviation
                          ess = F, 
                          rhat = F, 
                          conf.int = T,
                          conf.level = 0.95) %>% 
      #filter(str_detect(term, c('beta'))) %>% 
      dplyr::rename(mean=estimate,
                    std_deviation=std.error,
                    conf_low=conf.low,
                    conf_high=conf.high) %>% 
      mutate(cg = site_i,
             method_input_code = method_input_code_i,
             method_input_name = unique(sub_data$method_input_name),
             method = unique(sub_data$method),
             input_name = unique(sub_data$input_name),
             input_code = unique(sub_data$input_code),
             .before=1)     

}) %>% bind_rows()
 
}) %>% bind_rows()
  
coefftab %>% 
  saveRDS(file=here("outputs/ValidationCommonGarden/HeightModels/coefftab_M2.rds"))

4.2.3 Proportion of variance explained

The proportion of variance explained by the height models was estimated with the classical \(\mathcal{R}^{2}\), the residual-based Bayesian \(\mathcal{R}^{2}\) and the model-based Bayesian \(\mathcal{R}^{2}\).

The classical \(\mathcal{R}^{2}\) is calculated as follows:

\[\mathcal{R}^{2}=1 - \frac{Var(Y-\mu)}{Var(Y)} = \frac{Var(\mu)}{Var(Y)}\]

where \(\mu\) is the linear predictor of the model and \(Y\) the observed data (tree height).

The Bayesian version of \(\mathcal{R}^{2}\) introduced in is calculated as follows:

\[\mathcal{R}^{2}=Var_\mu/(Var_\mu + Var_ {res})\]

where \(Var_\mu\) is the variance of the modelled predictive means and \(Var_{res}\) is the residual variance.

In the residual-based Bayesian \(\mathcal{R}^{2}\), \(Var_{res}\) comes from draws from the residual distribution: \(Var_{res} = Var(Y-\mu)\). In the model-based Bayesian \(\mathcal{R}^{2}\), \(Var_{res}\) comes from the posterior quantities of the fitted model such as \(Var_{res} = \sigma^2\). As stated in , this Bayesian version of \(\mathcal{R}^{2}\) can be considered as a data-based estimate of the proportion of variance explained for new data.

More information on the Bayesian \(\mathcal{R}^{2}\) can also be found in the following vignette: Bayesian R2 and LOO-R2 by Aki Vehtari, Andrew Gelman, Ben Goodrich and Jonah Gabry.

In the main manuscript, we report results from the model-based Bayesian \(\mathcal{R}^{2}\).

4.2.4 Running model without predictor

Code
stancode = stan_model(here("scripts/StanModels/ValidationCommonGarden_HeightModel_M2_WithoutPredictor.stan"))

params_to_estimate <- c("beta_H","sigma","alpha_bloc","classic_R2","bayes_R2_mod","bayes_R2_res")

coefftab <- lapply(unique(height_data$cg), function(site_i){
  
    # Subset the data
    sub_data <- height_data %>% 
      filter(cg == site_i) %>% 
      left_join(pop_heights %>% dplyr::select(any_of(c("height", "pop"))) %>% dplyr::rename(init_height=height), by="pop")
      
    # Data in a list for Stan
    stanlist <- list(N = nrow(sub_data),
                     Y=(sub_data$height-mean(sub_data$height))/sd(sub_data$height),
                     H=(sub_data$init_height -mean(sub_data$init_height))/sd(sub_data$init_height),
                     nb_bloc = length(unique(sub_data$block)),
                     bloc = as.numeric(as.factor(sub_data$block)))

    # Run the models
    mod <- sampling(stancode, data = stanlist, iter = 2000, chains = 4, cores = 4, save_warmup = FALSE,
                    pars = params_to_estimate) 
    
    # Save the models
    mod %>% saveRDS(file=here(paste0("outputs/ValidationCommonGarden/HeightModels/stan_models/M2/M2_",site_i,"_WithoutPredictor.rds")))
    
    
    # Extract coefficients
    broom.mixed::tidyMCMC(mod,
                          pars=params_to_estimate,
                          droppars = NULL, 
                          robust = FALSE, # give mean and standard deviation
                          ess = F, 
                          rhat = F, 
                          conf.int = T,
                          conf.level = 0.95) %>% 
    dplyr::rename(mean=estimate,
                  std_deviation=std.error,
                  conf_low=conf.low,
                  conf_high=conf.high) %>% 
    mutate(cg=site_i,
           method="WithoutPredictor",
           .before=1)
    
    
  }) %>% bind_rows()

coefftab %>% saveRDS(file=here("outputs/ValidationCommonGarden/HeightModels/coefftab_M2_WithoutPredictor.rds"))

4.2.5 Model coefficients

We plot the mean and 95% credible intervals of the \(\beta_{X_1}\), \(\beta_{X_2}\) and \(\beta_H\) coefficients. Graph titles include the time in months corresponding to the age at which height and survival were recorded. Coefficients in the green area have the expected sign, reflecting lower height in populations with higher GO predictions.

Code
generate_interval_plots(model_i = "M2")
Joining with `by = join_by(cg, method_input_code, method_input_name, method, input_name, input_code)`
Warning in geom_rect(inherit.aes = FALSE, aes(xmin = -Inf, xmax = Inf, ymin = -Inf, : All aesthetics have length 1, but the data has 185 rows.
ℹ Please consider using `annotate()` or provide this layer with data containing a single row.
All aesthetics have length 1, but the data has 185 rows.
ℹ Please consider using `annotate()` or provide this layer with data containing a single row.
Warning in geom_rect(inherit.aes = FALSE, aes(xmin = -Inf, xmax = Inf, ymin = -Inf, : All aesthetics have length 1, but the data has 150 rows.
ℹ Please consider using `annotate()` or provide this layer with data containing a single row.
Joining with `by = join_by(cg, method_input_code, method_input_name, method, input_name, input_code)`
Joining with `by = join_by(cg, method_input_code, method_input_name, method, input_name, input_code)`
[[1]]
Warning in geom_rect(inherit.aes = FALSE, aes(xmin = -Inf, xmax = Inf, ymin = -Inf, : All aesthetics have length 1, but the data has 185 rows.
ℹ Please consider using `annotate()` or provide this layer with data containing a single row.


[[2]]


[[3]]


[[4]]


[[5]]


[[6]]

4.2.6 Vizualising predicted associations

4.2.6.1 Mean height predictions

We plot the mean predicted relationships between tree height and GO predictions / CTD.

Code
generate_poly_plots(model_i = "M2")

4.2.6.2 Height predictions with uncertainty intervals

The graphs above do not show the estimate uncertainty. We generate some graphs to visualize the uncertainty around the estimates :

  • the uncertainty in the mean estimates (i.e. variability in posterior draws of the linear predictor \(\mu\))

  • the uncertainty in tree height predictions (i.e. after accounting for \(\sigma\) in the predictions).

Code
#############################################################################
# Visualizing the relationships between GO predictions and height predictions
#############################################################################
model_i <- "M2"

coefftab <- readRDS(file=here(paste0("outputs/ValidationCommonGarden/HeightModels/coefftab_",model_i,".rds"))) %>% 
  mutate(cg_name = case_when(cg == "asturias" ~ "Asturias, Spain (37 months)",
                             cg == "bordeaux" ~ "Bordeaux, France (85 months)",
                             cg == "caceres" ~ "Cáceres, Spain (8 months)",
                             cg == "madrid" ~ "Madrid, Spain (13 months)",
                             cg == "portugal" ~ "Fundão, Portugal (27 months)"))

lapply(unique(coefftab$cg), function(cg_i){
  
  lapply(unique(coefftab$method_input_code), function(method_input_code_i){
    
    sub_coefftab <- coefftab %>% 
      filter(method_input_code == method_input_code_i & cg == cg_i)
    
    # Load the stanlist
    stanlist <- readRDS(file = here(paste0("outputs/ValidationCommonGarden/HeightModels/stan_models/M2/",cg_i,"_",method_input_code_i,".rds")))[[2]]
  
    # Loading the model
    mod <- readRDS(file = here(paste0("outputs/ValidationCommonGarden/HeightModels/stan_models/M2/",cg_i,"_",method_input_code_i,".rds")))[[1]]

    # we extract the posterior distributions of all parameters
    post <- as.data.frame(mod)

    # Vector of predictor values (based on the min and max of the GO predictions/CTD )
    x_vec <- seq(min(stanlist$X),max(stanlist$X),0.05)

    # Predicting the linear predictor mu (predicting mean-centered height without sigma uncertainty)
    ################################################################################################

    # we extract the posterior draws of mu, and its mean and HDPIs for each predictor value

    # Function to predict mean-centered height in block 2 and with the initial tree height fixed to the mean (i.e. fixed to zero as explanatory variables were mean-centered)
    f_mu <- function(x) post$`alpha_bloc[2]` + post$`beta_X1` * x + post$`beta_X2` * x^2 + post$`beta_H` * mean(stanlist$H)#mean(sub_data$init_height) 

    mu_pred <- 
      sapply(x_vec, f_mu) %>% 
      tibble::as_tibble() %>%
      rename_all(function(x) x_vec) %>%
      mutate(Iter = row_number()) %>%
      gather(x, y, -Iter) %>%
      group_by(x) %>%
      mutate(hpdi_l = HDInterval::hdi(y, credMass = 0.90)[1],
             hpdi_h = HDInterval::hdi(y, credMass = 0.90)[2]) %>%
      mutate(mu_mean = mean(y)) %>%
      ungroup() %>%
      mutate(x = as.numeric(x))
    
    # Keep mean and 90% HDPI of mu (one value for each iteration)
    mu_mean_df <- mu_pred %>%
      dplyr::select(x,mu_mean,hpdi_l,hpdi_h) %>% 
      dplyr::distinct()

   # Predicting mean-centered height with sigma uncertainty
   #########################################################

   # we extract the posterior draws of height predictions, and its mean and PIs for each predictor value
    
    y_pred <- 
      sapply(x_vec,
         function(x)
           rnorm(NROW(post),
                 post$`alpha_bloc[2]` + post$`beta_X1` * x + post$`beta_X2` * x^2 + post$`beta_H` *  mean(stanlist$H),#mean(sub_data$init_height),
                 post$sigma)) %>%
  as_tibble() %>%
  rename_all(function(x) x_vec) %>%
  mutate(Iter = row_number()) %>%
  gather(x, y, -Iter) %>%
  group_by(x) %>%
  mutate(hpdi_l = HDInterval::hdi(y, credMass = 0.90)[1],
         hpdi_h = HDInterval::hdi(y, credMass = 0.90)[2]
         # pi_l = rethinking::PI(y, prob = 0.90)[1],
         # pi_h = rethinking::PI(y, prob = 0.90)[2]
         ) %>%
  ungroup() %>%
  mutate(x = as.numeric(x)) %>% 
  dplyr::select(x,hpdi_l,hpdi_h) %>% 
  dplyr::distinct()


  # Plots
  #######

  # Plot options 
  y_limits <- c(-3,3)
  if(unique(sub_coefftab$method == "CTD")){x_axis <- "Mean-centered CTD"} else {
    x_axis <- "Mean-centered GO predictions"}
  y_axis <- "Mean-centered tree height"

  p <- ggplot() +
    ylim(y_limits) +
    labs(y=y_axis, x=x_axis) +
    theme_bw()

  # First 100 posterior draws of the linear predictor mu
  p1 <- p +
    geom_point(data = mu_pred %>% filter(Iter < 101),
               aes(x, y), alpha = .2, color = 'dodgerblue')
  
  # Mean (line) and 90% HPDI (shaded region) of the linear predictor mu
  p2 <- p +
    geom_ribbon(data = y_pred,
                mapping = aes(x=x, ymin=hpdi_l, ymax=hpdi_h), alpha = .1, fill='dodgerblue') +
    geom_ribbon(data = mu_mean_df,
                aes(x = x, ymin = hpdi_l, ymax = hpdi_h),
                alpha = .2,
                fill='dodgerblue') +
    geom_line(data = mu_mean_df,
              aes(x=x, y=mu_mean), 
              col="dodgerblue4")

  p1_p2 <- ggarrange(p1, p2 + labs(y=""), nrow = 1)

  # Title
  title <- ggdraw() + 
    draw_label(paste0(unique(sub_coefftab$cg_name)," - ",unique(sub_coefftab$method_input_name)),
               fontface = 'bold',
               x = 0, 
               hjust = 0) +   
  theme(plot.margin = margin(0, 0, 0, 7))

  # merge title and plots
  p1_p2 <- plot_grid(title, p1_p2, ncol = 1, rel_heights = c(0.1, 1))

  ggsave(p1_p2,
         file=here(paste0("figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/",model_i,
                        "_",cg_i,"_",method_input_code_i,".pdf")),
         device="pdf",
         height=6,
         width=11)
  })
  })
[[1]]
[[1]][[1]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_asturias_CTD_bio1.pdf"

[[1]][[2]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_asturias_CTD_bio12.pdf"

[[1]][[3]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_asturias_CTD_bio15.pdf"

[[1]][[4]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_asturias_CTD_bio3.pdf"

[[1]][[5]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_asturias_CTD_bio4.pdf"

[[1]][[6]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_asturias_CTD_SHM.pdf"

[[1]][[7]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_asturias_CTD_EucliDist.pdf"

[[1]][[8]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_asturias_GDM_all_cand.pdf"

[[1]][[9]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_asturias_GF_all_cand.pdf"

[[1]][[10]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_asturias_LFMM_all_cand.pdf"

[[1]][[11]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_asturias_pRDA_all_cand.pdf"

[[1]][[12]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_asturias_RDA_all_cand.pdf"

[[1]][[13]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_asturias_GDM_common_cand.pdf"

[[1]][[14]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_asturias_GF_common_cand.pdf"

[[1]][[15]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_asturias_LFMM_common_cand.pdf"

[[1]][[16]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_asturias_pRDA_common_cand.pdf"

[[1]][[17]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_asturias_RDA_common_cand.pdf"

[[1]][[18]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_asturias_GDM_corrected_cand.pdf"

[[1]][[19]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_asturias_GF_corrected_cand.pdf"

[[1]][[20]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_asturias_LFMM_corrected_cand.pdf"

[[1]][[21]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_asturias_pRDA_corrected_cand.pdf"

[[1]][[22]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_asturias_RDA_corrected_cand.pdf"

[[1]][[23]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_asturias_GDM_randomfreq_control.pdf"

[[1]][[24]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_asturias_GF_randomfreq_control.pdf"

[[1]][[25]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_asturias_LFMM_randomfreq_control.pdf"

[[1]][[26]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_asturias_pRDA_randomfreq_control.pdf"

[[1]][[27]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_asturias_RDA_randomfreq_control.pdf"

[[1]][[28]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_asturias_GDM_simfreq_control.pdf"

[[1]][[29]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_asturias_GF_simfreq_control.pdf"

[[1]][[30]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_asturias_LFMM_simfreq_control.pdf"

[[1]][[31]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_asturias_pRDA_simfreq_control.pdf"

[[1]][[32]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_asturias_RDA_simfreq_control.pdf"

[[1]][[33]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_asturias_GDM_all.pdf"

[[1]][[34]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_asturias_GF_all.pdf"

[[1]][[35]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_asturias_LFMM_all.pdf"

[[1]][[36]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_asturias_pRDA_all.pdf"

[[1]][[37]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_asturias_RDA_all.pdf"


[[2]]
[[2]][[1]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_bordeaux_CTD_bio1.pdf"

[[2]][[2]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_bordeaux_CTD_bio12.pdf"

[[2]][[3]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_bordeaux_CTD_bio15.pdf"

[[2]][[4]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_bordeaux_CTD_bio3.pdf"

[[2]][[5]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_bordeaux_CTD_bio4.pdf"

[[2]][[6]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_bordeaux_CTD_SHM.pdf"

[[2]][[7]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_bordeaux_CTD_EucliDist.pdf"

[[2]][[8]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_bordeaux_GDM_all_cand.pdf"

[[2]][[9]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_bordeaux_GF_all_cand.pdf"

[[2]][[10]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_bordeaux_LFMM_all_cand.pdf"

[[2]][[11]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_bordeaux_pRDA_all_cand.pdf"

[[2]][[12]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_bordeaux_RDA_all_cand.pdf"

[[2]][[13]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_bordeaux_GDM_common_cand.pdf"

[[2]][[14]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_bordeaux_GF_common_cand.pdf"

[[2]][[15]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_bordeaux_LFMM_common_cand.pdf"

[[2]][[16]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_bordeaux_pRDA_common_cand.pdf"

[[2]][[17]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_bordeaux_RDA_common_cand.pdf"

[[2]][[18]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_bordeaux_GDM_corrected_cand.pdf"

[[2]][[19]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_bordeaux_GF_corrected_cand.pdf"

[[2]][[20]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_bordeaux_LFMM_corrected_cand.pdf"

[[2]][[21]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_bordeaux_pRDA_corrected_cand.pdf"

[[2]][[22]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_bordeaux_RDA_corrected_cand.pdf"

[[2]][[23]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_bordeaux_GDM_randomfreq_control.pdf"

[[2]][[24]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_bordeaux_GF_randomfreq_control.pdf"

[[2]][[25]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_bordeaux_LFMM_randomfreq_control.pdf"

[[2]][[26]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_bordeaux_pRDA_randomfreq_control.pdf"

[[2]][[27]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_bordeaux_RDA_randomfreq_control.pdf"

[[2]][[28]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_bordeaux_GDM_simfreq_control.pdf"

[[2]][[29]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_bordeaux_GF_simfreq_control.pdf"

[[2]][[30]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_bordeaux_LFMM_simfreq_control.pdf"

[[2]][[31]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_bordeaux_pRDA_simfreq_control.pdf"

[[2]][[32]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_bordeaux_RDA_simfreq_control.pdf"

[[2]][[33]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_bordeaux_GDM_all.pdf"

[[2]][[34]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_bordeaux_GF_all.pdf"

[[2]][[35]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_bordeaux_LFMM_all.pdf"

[[2]][[36]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_bordeaux_pRDA_all.pdf"

[[2]][[37]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_bordeaux_RDA_all.pdf"


[[3]]
[[3]][[1]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_caceres_CTD_bio1.pdf"

[[3]][[2]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_caceres_CTD_bio12.pdf"

[[3]][[3]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_caceres_CTD_bio15.pdf"

[[3]][[4]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_caceres_CTD_bio3.pdf"

[[3]][[5]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_caceres_CTD_bio4.pdf"

[[3]][[6]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_caceres_CTD_SHM.pdf"

[[3]][[7]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_caceres_CTD_EucliDist.pdf"

[[3]][[8]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_caceres_GDM_all_cand.pdf"

[[3]][[9]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_caceres_GF_all_cand.pdf"

[[3]][[10]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_caceres_LFMM_all_cand.pdf"

[[3]][[11]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_caceres_pRDA_all_cand.pdf"

[[3]][[12]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_caceres_RDA_all_cand.pdf"

[[3]][[13]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_caceres_GDM_common_cand.pdf"

[[3]][[14]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_caceres_GF_common_cand.pdf"

[[3]][[15]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_caceres_LFMM_common_cand.pdf"

[[3]][[16]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_caceres_pRDA_common_cand.pdf"

[[3]][[17]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_caceres_RDA_common_cand.pdf"

[[3]][[18]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_caceres_GDM_corrected_cand.pdf"

[[3]][[19]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_caceres_GF_corrected_cand.pdf"

[[3]][[20]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_caceres_LFMM_corrected_cand.pdf"

[[3]][[21]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_caceres_pRDA_corrected_cand.pdf"

[[3]][[22]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_caceres_RDA_corrected_cand.pdf"

[[3]][[23]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_caceres_GDM_randomfreq_control.pdf"

[[3]][[24]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_caceres_GF_randomfreq_control.pdf"

[[3]][[25]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_caceres_LFMM_randomfreq_control.pdf"

[[3]][[26]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_caceres_pRDA_randomfreq_control.pdf"

[[3]][[27]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_caceres_RDA_randomfreq_control.pdf"

[[3]][[28]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_caceres_GDM_simfreq_control.pdf"

[[3]][[29]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_caceres_GF_simfreq_control.pdf"

[[3]][[30]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_caceres_LFMM_simfreq_control.pdf"

[[3]][[31]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_caceres_pRDA_simfreq_control.pdf"

[[3]][[32]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_caceres_RDA_simfreq_control.pdf"

[[3]][[33]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_caceres_GDM_all.pdf"

[[3]][[34]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_caceres_GF_all.pdf"

[[3]][[35]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_caceres_LFMM_all.pdf"

[[3]][[36]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_caceres_pRDA_all.pdf"

[[3]][[37]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_caceres_RDA_all.pdf"


[[4]]
[[4]][[1]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_madrid_CTD_bio1.pdf"

[[4]][[2]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_madrid_CTD_bio12.pdf"

[[4]][[3]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_madrid_CTD_bio15.pdf"

[[4]][[4]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_madrid_CTD_bio3.pdf"

[[4]][[5]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_madrid_CTD_bio4.pdf"

[[4]][[6]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_madrid_CTD_SHM.pdf"

[[4]][[7]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_madrid_CTD_EucliDist.pdf"

[[4]][[8]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_madrid_GDM_all_cand.pdf"

[[4]][[9]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_madrid_GF_all_cand.pdf"

[[4]][[10]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_madrid_LFMM_all_cand.pdf"

[[4]][[11]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_madrid_pRDA_all_cand.pdf"

[[4]][[12]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_madrid_RDA_all_cand.pdf"

[[4]][[13]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_madrid_GDM_common_cand.pdf"

[[4]][[14]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_madrid_GF_common_cand.pdf"

[[4]][[15]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_madrid_LFMM_common_cand.pdf"

[[4]][[16]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_madrid_pRDA_common_cand.pdf"

[[4]][[17]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_madrid_RDA_common_cand.pdf"

[[4]][[18]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_madrid_GDM_corrected_cand.pdf"

[[4]][[19]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_madrid_GF_corrected_cand.pdf"

[[4]][[20]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_madrid_LFMM_corrected_cand.pdf"

[[4]][[21]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_madrid_pRDA_corrected_cand.pdf"

[[4]][[22]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_madrid_RDA_corrected_cand.pdf"

[[4]][[23]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_madrid_GDM_randomfreq_control.pdf"

[[4]][[24]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_madrid_GF_randomfreq_control.pdf"

[[4]][[25]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_madrid_LFMM_randomfreq_control.pdf"

[[4]][[26]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_madrid_pRDA_randomfreq_control.pdf"

[[4]][[27]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_madrid_RDA_randomfreq_control.pdf"

[[4]][[28]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_madrid_GDM_simfreq_control.pdf"

[[4]][[29]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_madrid_GF_simfreq_control.pdf"

[[4]][[30]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_madrid_LFMM_simfreq_control.pdf"

[[4]][[31]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_madrid_pRDA_simfreq_control.pdf"

[[4]][[32]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_madrid_RDA_simfreq_control.pdf"

[[4]][[33]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_madrid_GDM_all.pdf"

[[4]][[34]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_madrid_GF_all.pdf"

[[4]][[35]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_madrid_LFMM_all.pdf"

[[4]][[36]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_madrid_pRDA_all.pdf"

[[4]][[37]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_madrid_RDA_all.pdf"


[[5]]
[[5]][[1]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_portugal_CTD_bio1.pdf"

[[5]][[2]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_portugal_CTD_bio12.pdf"

[[5]][[3]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_portugal_CTD_bio15.pdf"

[[5]][[4]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_portugal_CTD_bio3.pdf"

[[5]][[5]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_portugal_CTD_bio4.pdf"

[[5]][[6]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_portugal_CTD_SHM.pdf"

[[5]][[7]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_portugal_CTD_EucliDist.pdf"

[[5]][[8]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_portugal_GDM_all_cand.pdf"

[[5]][[9]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_portugal_GF_all_cand.pdf"

[[5]][[10]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_portugal_LFMM_all_cand.pdf"

[[5]][[11]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_portugal_pRDA_all_cand.pdf"

[[5]][[12]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_portugal_RDA_all_cand.pdf"

[[5]][[13]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_portugal_GDM_common_cand.pdf"

[[5]][[14]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_portugal_GF_common_cand.pdf"

[[5]][[15]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_portugal_LFMM_common_cand.pdf"

[[5]][[16]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_portugal_pRDA_common_cand.pdf"

[[5]][[17]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_portugal_RDA_common_cand.pdf"

[[5]][[18]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_portugal_GDM_corrected_cand.pdf"

[[5]][[19]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_portugal_GF_corrected_cand.pdf"

[[5]][[20]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_portugal_LFMM_corrected_cand.pdf"

[[5]][[21]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_portugal_pRDA_corrected_cand.pdf"

[[5]][[22]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_portugal_RDA_corrected_cand.pdf"

[[5]][[23]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_portugal_GDM_randomfreq_control.pdf"

[[5]][[24]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_portugal_GF_randomfreq_control.pdf"

[[5]][[25]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_portugal_LFMM_randomfreq_control.pdf"

[[5]][[26]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_portugal_pRDA_randomfreq_control.pdf"

[[5]][[27]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_portugal_RDA_randomfreq_control.pdf"

[[5]][[28]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_portugal_GDM_simfreq_control.pdf"

[[5]][[29]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_portugal_GF_simfreq_control.pdf"

[[5]][[30]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_portugal_LFMM_simfreq_control.pdf"

[[5]][[31]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_portugal_pRDA_simfreq_control.pdf"

[[5]][[32]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_portugal_RDA_simfreq_control.pdf"

[[5]][[33]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_portugal_GDM_all.pdf"

[[5]][[34]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_portugal_GF_all.pdf"

[[5]][[35]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_portugal_LFMM_all.pdf"

[[5]][[36]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_portugal_pRDA_all.pdf"

[[5]][[37]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/HeightPredictions_MuPosteriorDraws_UncertaintyIntervals/M2_portugal_RDA_all.pdf"

4.2.7 Correlation coefficients

Stan code of the height model without the genomic offset predictions or the CTD (i.e. only with height as predictor).

Code
stancode = stan_model(here("scripts/StanModels/ValidationCommonGarden_HeightModel_M2_WithoutPredictor_WithPredictions.stan"))
print(stancode)
Code
lapply(unique(height_data$cg), function(site_i){
  
  # Subset the data
  sub_data <- height_data %>% 
    filter(cg == site_i) %>% 
    left_join(pop_heights %>% dplyr::select(any_of(c("height", "pop"))) %>% dplyr::rename(init_height=height), by="pop")
  
  # Data in a list for Stan
  stanlist <- list(N = nrow(sub_data),
                   Y=(sub_data$height-mean(sub_data$height))/sd(sub_data$height),
                   H=(sub_data$init_height -mean(sub_data$init_height))/sd(sub_data$init_height),
                   nb_bloc = length(unique(sub_data$block)),
                   bloc = as.numeric(as.factor(sub_data$block)))

  # Run the models
  mod <- sampling(stancode, 
                  data = stanlist, 
                  iter = 2000, 
                  chains = 4, 
                  cores = 4, 
                  save_warmup = FALSE,
                  warmup = 1400,
                  pars = "y_pred")
  
  # Save the stanlist and the model
  saveRDS(list(data=sub_data, mod=mod),
          file=here(paste0("outputs/ValidationCommonGarden/HeightModels/stan_models/M2/M2_WithoutPredictorWithPredictions_",site_i,".rds")))
                  
})
Code
corrtab <-lapply(unique(height_data$cg), function(site_i){
  
  mod <- readRDS(mod, file=here(paste0("outputs/ValidationCommonGarden/HeightModels/stan_models/M2/M2_WithoutPredictorWithPredictions_",site_i,".rds")))
  
  # Model coefficients
  coefftab <- broom.mixed::tidyMCMC(mod$mod,
                                    pars= "y_pred",
                                    droppars = NULL, 
                                    robust = FALSE, # give mean and standard deviation
                                    ess = T, # effective sample size estimates
                                    rhat = T, # Rhat estimates
                                    conf.int = T, 
                                    conf.level = 0.95 # include 95% credible intervals
                                    ) %>% 
    dplyr::rename(conf95_low=conf.low,
                  conf95_high=conf.high,
                  mean=estimate,
                  std_deviation=std.error)  
  
  res_height_df <- mod$data %>%
    mutate(height_sd = (height-mean(height))/sd(height)) %>% 
    bind_cols(coefftab) %>% 
    mutate(res_height = height_sd - mean) %>% # observed height - predicted height 
    dplyr::select(cg, pop, res_height) %>% 
    group_by(pop) %>% 
    summarise(mean_res_height_pop =mean(res_height))
  
  
  
  lapply(unique(df$method_input_code), function(method_input_code_i){
    
    sub_data <- df %>% 
      filter(method_input_code == method_input_code_i & cg == site_i) %>% 
      inner_join(res_height_df, by="pop")
    
    tibble(pearson_correlation = cor.test(sub_data$varX, sub_data$mean_res_height_pop, method="pearson")$estimate[[1]],
           pearson_pvalue = cor.test(sub_data$varX, sub_data$mean_res_height_pop, method="pearson")$p.value,
           spearman_correlation = cor.test(sub_data$varX, sub_data$mean_res_height_pop, method="spearman")$estimate[[1]],
           spearman_pvalue = cor.test(sub_data$varX, sub_data$mean_res_height_pop, method="spearman")$p.value,
           method_input_code = method_input_code_i,
           method_input_name = unique(sub_data$method_input_name),
           input_name = unique(sub_data$input_name),
           input_code = unique(sub_data$input_code),
           method = unique(sub_data$method),
           cg = site_i)
    }) %>% bind_rows()   
  }) %>% bind_rows()

saveRDS(corrtab, here("outputs/ValidationCommonGarden/HeightModels/corrtab_M2.rds"))
Code
corrtab <- readRDS(here("outputs/ValidationCommonGarden/HeightModels/corrtab_M2.rds"))

correlation_types <- c("pearson_correlation", "spearman_correlation")

lapply(correlation_types, function(coeff){
  
  p <- corrtab %>%
    pivot_longer(values_to = "correlation_estimate", names_to = "correlation_type", cols = all_of(correlation_types)) %>% 
    mutate(correlation_pvalue = ifelse(correlation_type == "pearson_correlation", pearson_pvalue, spearman_pvalue)) %>% 
    dplyr::select(-pearson_pvalue,-spearman_pvalue) %>% 
    mutate(pvalue_signi = ifelse(correlation_pvalue < 0.05, "p-value < 0.05", "p-value ≥ 0.05")) %>% 
    filter(correlation_type == coeff) %>% 
    mutate(cg_name = case_when(cg == "asturias" ~ "Asturias, Spain (37 months)",
                               cg == "bordeaux" ~ "Bordeaux, France (85 months)",
                               cg == "caceres" ~ "Cáceres, Spain (8 months)",
                               cg == "madrid" ~ "Madrid, Spain (13 months)",
                               cg == "portugal" ~ "Fundão, Portugal (27 months)")) %>% 
    mutate(input_name = factor(input_name, levels = c(unique(corrtab$input_name)[[13]],
                                                      unique(corrtab$input_name)[[11]],
                                                      unique(corrtab$input_name)[[12]],
                                                      unique(corrtab$input_name)[[10]],
                                                      unique(corrtab$input_name)[[8]],
                                                      unique(corrtab$input_name)[[9]],
                                                      unique(corrtab$input_name)[[1]],
                                                      unique(corrtab$input_name)[[2]],
                                                      unique(corrtab$input_name)[[3]],
                                                      unique(corrtab$input_name)[[4]],
                                                      unique(corrtab$input_name)[[5]],
                                                      unique(corrtab$input_name)[[6]],
                                                      unique(corrtab$input_name)[[7]])),
           method = factor(method, levels = c(unique(corrtab$method)[[2]],
                                              unique(corrtab$method)[[3]],
                                              unique(corrtab$method)[[4]],
                                              unique(corrtab$method)[[5]],
                                              unique(corrtab$method)[[6]],
                                              unique(corrtab$method)[[1]])),
           cg_name = factor(cg_name, levels = c("Madrid, Spain (13 months)",
                                                "Cáceres, Spain (8 months)",
                                                "Asturias, Spain (37 months)",
                                                "Bordeaux, France (85 months)",
                                                "Fundão, Portugal (27 months)")))

# Plots with CTD and SNP sets
# ===========================
p_allvar <- p %>% ggplot(aes(x = method,
                             y = correlation_estimate,
                             color = input_name,
                             shape = pvalue_signi)) +
  geom_rect(inherit.aes=FALSE, aes(xmin=-Inf, xmax=Inf, ymin=-Inf, ymax=0), color="transparent", fill="green", alpha=0.005) + 
  geom_hline(yintercept = 0,color="gray") +
  geom_point(position = position_dodge(width = .5), size=5) +
  facet_wrap(~cg_name, ncol=2) + 
  {if(coeff=="pearson_correlation") ylab("Pearson correlation estimates")} + 
  {if(coeff=="spearman_correlation") ylab("Spearman correlation estimates (rank correlation)") } + 
  xlab("") +
  #ylim(c(-0.44,0.55)) +
  scale_color_manual(values=colors_coeff)+
  scale_shape_manual(values = c(16,18)) + #, name="p-value of the correlations"
  theme_bw() +
  labs(color="",shape="") +
  theme(axis.text.x =  element_text(size=13),
        axis.text.y = element_text(size=13),
        axis.title.y = element_text(size=16),
        legend.title = element_blank(), 
        strip.text.x = element_text(size = 16),
        strip.background = element_blank(),
        legend.position = c(0.77,0.15),
        legend.text=element_text(size=13),
        panel.grid.minor.x=element_blank(),
        panel.grid.major.x=element_blank()) +
  guides(color = guide_legend(ncol=1),
         shape = guide_legend(override.aes = list(size =2 )))
p_allvar
# save in pdf and png
ggsave(p_allvar,
       file=here(paste0("figs/ValidationCommonGarden/HeightModels/Correlations/",coeff,"_SNPsandCTD.pdf")),
       device="pdf",
       height=13.6,
       width=15)

ggsave(p_allvar,
       file=here(paste0("figs/ValidationCommonGarden/HeightModels/Correlations/",coeff,"_SNPsandCTD.png")),
       height=13.6,
       width=15)



# Adding some images to represent the climatic differences among common gardens
annotation_custom2 <- function (grob, xmin = -Inf, xmax = Inf, ymin = -Inf, ymax = Inf, data){ layer(data = data, stat = StatIdentity, position = PositionIdentity, 
        geom = ggplot2:::GeomCustomAnn,
        inherit.aes = TRUE, params = list(grob = grob, 
                                          xmin = xmin, xmax = xmax, 
                                          ymin = ymin, ymax = ymax))}

df_png <- p %>% 
  dplyr::select(cg_name) %>%
  distinct %>% 
  dplyr::mutate(image = case_when(cg_name == "Asturias, Spain (37 months)" ~ "reports/cloud.png",
                                  cg_name == "Bordeaux, France (85 months)" ~ "reports/cloud.png",
                                  cg_name == "Cáceres, Spain (8 months)" ~ "reports/sun.png",
                                  cg_name == "Madrid, Spain (13 months)" ~ "reports/sun.png",
                                  cg_name == "Fundão, Portugal (27 months)" ~ 'reports/cloud_and_sun.png'))

list_pngs <- lapply(unique(p$cg_name), function(cg_name_i){
  
sub <-  p %>% 
  filter(cg_name==cg_name_i) %>%
  slice(1)

png_cg = annotation_custom2(rasterGrob(readPNG(here(df_png[df_png$cg_name==cg_name_i,"image"])),interpolate=TRUE), 
                            ymin = -0.65,
                            ymax= -0.45,
                            xmin = 0.5,
                            xmax = 1.1, 
                            data=sub)
})

p_allvar_images <- p_allvar + list_pngs[[1]]+ list_pngs[[2]]+ list_pngs[[3]]+ list_pngs[[4]]+ list_pngs[[5]]

p_allvar_images

# save in pdf and png
ggsave(p_allvar_images,
       file=here(paste0("figs/ValidationCommonGarden/HeightModels/Correlations/",coeff,"_SNPsandCTD_WithImages.pdf")),
       device="pdf",
       height=13.6,
       width=15)

})
Warning in geom_rect(inherit.aes = FALSE, aes(xmin = -Inf, xmax = Inf, ymin = -Inf, : All aesthetics have length 1, but the data has 185 rows.
ℹ Please consider using `annotate()` or provide this layer with data containing a single row.
Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, : for 'p-value ≥ 0.05' in 'mbcsToSbcs': >= substituted for ≥ (U+2265)
Warning in geom_rect(inherit.aes = FALSE, aes(xmin = -Inf, xmax = Inf, ymin = -Inf, : All aesthetics have length 1, but the data has 185 rows.
ℹ Please consider using `annotate()` or provide this layer with data containing a single row.
All aesthetics have length 1, but the data has 185 rows.
ℹ Please consider using `annotate()` or provide this layer with data containing a single row.
Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, : for 'p-value ≥ 0.05' in 'mbcsToSbcs': >= substituted for ≥ (U+2265)
Warning in geom_rect(inherit.aes = FALSE, aes(xmin = -Inf, xmax = Inf, ymin = -Inf, : All aesthetics have length 1, but the data has 185 rows.
ℹ Please consider using `annotate()` or provide this layer with data containing a single row.
Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, : for 'p-value ≥ 0.05' in 'mbcsToSbcs': >= substituted for ≥ (U+2265)
Warning in geom_rect(inherit.aes = FALSE, aes(xmin = -Inf, xmax = Inf, ymin = -Inf, : All aesthetics have length 1, but the data has 185 rows.
ℹ Please consider using `annotate()` or provide this layer with data containing a single row.
All aesthetics have length 1, but the data has 185 rows.
ℹ Please consider using `annotate()` or provide this layer with data containing a single row.
Warning in grid.Call.graphics(C_text, as.graphicsAnnot(x$label), x$x, x$y, : for 'p-value ≥ 0.05' in 'mbcsToSbcs': >= substituted for ≥ (U+2265)
[[1]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/Correlations/pearson_correlation_SNPsandCTD_WithImages.pdf"

[[2]]
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/figs/ValidationCommonGarden/HeightModels/Correlations/spearman_correlation_SNPsandCTD_WithImages.pdf"

4.3 Model 3

4.3.1 Model equation and code

Model 3 = Model 2 + Initial height as a confounder + Clone fixed intercepts

\[\begin{align*} Y_{ipb} &\sim \mathcal{N}(\mu_{pb},\sigma^{2}) \\ \mu_{pb} &= \beta_0 + B_b + \beta_{X1}X_p + \beta_{X2}X^2_p + \beta_H H_p\\ \beta_0 &\sim \mathcal{N}(\text{mean}(Y),2) \\ \sigma &\sim \text{Exp}(1) \\ \begin{bmatrix} B_b \\ C_c \\ \beta_{X1} \\ \beta_{X2} \\ \beta_H \end{bmatrix} &\sim \mathcal{N}(0,1) \\ \end{align*}\]

  • \(Y_{ipb}\) is the height of the individual \(i\) in the population \(p\) and block \(b\).
  • \(\sigma^{2}\) is the residual variance of the model.
  • \(B_b\) are the block intercepts.
  • \(X_p\) is the GO or CTD of the population \(p\), with \(\beta_{X1}\) and \(\beta_{X2}\) being its linear and quadratic coefficients, respectively. The quadratic term was included to allow for potential nonlinearity in the response, following Fitzpatrick et al. (2021).
  • \(H_p\) is a proxy of the initial height of the trees from population \(p\), i.e. when trees were planted. For that, we used the population varying intercepts calculated across all common gardens in the model 1 of Archambeau et al. (2022).
  • \(\beta_0\) is the model global intercept.
  • \(C_c\) are the clone intercepts.
Code
stancode = stan_model(here("scripts/StanModels/ValidationCommonGarden_HeightModel_M3.stan"))
print(stancode)
S4 class stanmodel 'anon_model' coded as follows:
data {
  int N;
  vector[N] Y;                              // Response variable (individual height)
  vector[N] X;                              // Genomic offset or climatic transfer distances
  vector[N] H;                              // Initial height of the populations
  int<lower=0> nb_bloc;                     // Number of blocks
  int<lower=0, upper=nb_bloc> bloc[N];      // Blocks
  int<lower=0> nb_clon;                     // Number of clones
  int<lower=0, upper=nb_clon> clon[N];      // Clones
}
parameters {
  vector[nb_bloc] alpha_bloc;               // Intercepts of the blocks
  vector[nb_clon] alpha_clon;               // Intercepts of the clones
  real beta_0;                              // Global intercept
  real beta_X1;                             // Linear coefficent of GO or CTD
  real beta_X2;                             // Quadratic coefficent of GO or CTD
  real beta_H;                              // Linear coefficient of the initial height of the populations
  real<lower = 0>  sigma;                   // Residual variance of the model
}
transformed parameters {
  vector[N] mu;                             // Linear predictor
  real R_squared;                           // R^2 to evaluate the goodness of fit of the model

  mu = beta_0 + alpha_bloc[bloc] + alpha_clon[clon] + beta_X1 * X + beta_X2 * square(X) + beta_H * H;
  R_squared = 1 - variance(Y - mu) / variance(Y);
}
model {

  // Likelihood
  Y ~ normal(mu, sigma);

  // Priors
  sigma ~ exponential(1);
  alpha_bloc ~ std_normal();
  alpha_clon ~ std_normal();
  beta_0 ~ std_normal();
  beta_X1 ~ std_normal();
  beta_X2 ~ std_normal();
  beta_H ~ std_normal();
}
// generated quantities{
//   // log likelihood for loo
//   vector[N] log_lik;
//   vector[N] muhat;
//   for (n in 1:N) {
//     log_lik[n] = normal_lpdf(Y[n] |mu[n],sigma);
//     muhat[n] = normal_rng(mu[n], sigma);
//   }
// } 
Code
params_to_estimate <- c("beta_X1","beta_X2","beta_H","R_squared","sigma")

4.3.2 Run the models

Code
lapply(unique(height_data$cg), function(site_i){
  
  lapply(unique(df$method_input_code), function(method_input_code_i){
  
    # Subset the data
    df_sub <- df %>% filter(method_input_code == method_input_code_i & cg == site_i)
    
    sub_data <- height_data %>% 
      filter(cg == site_i) %>% 
      left_join(df_sub, by = c("cg","pop")) %>% 
      left_join(pop_heights %>% dplyr::select(any_of(c("height", "pop"))) %>% dplyr::rename(init_height=height), by="pop")
      
    # Data in a list for Stan
    stanlist <- list(N = nrow(sub_data),
                     Y=(sub_data$height-mean(sub_data$height))/sd(sub_data$height),
                     X=(sub_data$varX -mean(sub_data$varX))/sd(sub_data$varX),
                     H=(sub_data$init_height -mean(sub_data$init_height))/sd(sub_data$init_height),
                     nb_bloc = length(unique(sub_data$block)),
                     nb_clon = length(unique(sub_data$clon)),
                     bloc = as.numeric(as.factor(sub_data$block)),
                     clon = as.numeric(as.factor(sub_data$clon)))

    # Run the models
    mod <- sampling(stancode, 
                    data = stanlist, 
                    iter = 2000, 
                    warmup = 1400, 
                    chains = 4, 
                    cores = 4, 
                    save_warmup = FALSE,
                    pars = params_to_estimate) 
    
    # Save the model and the stanlist
    list(mod = mod, 
         stanlist = stanlist) %>% 
      saveRDS(file=here(paste0("outputs/ValidationCommonGarden/HeightModels/stan_models/M3/",site_i,"_",method_input_code_i,".rds")))
    
    
  })
})

Warning message: “Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable. Running the chains for more iterations may help. See https://mc-stan.org/misc/warnings.html#bulk-ess”

Code
coefftab <- lapply(unique(height_data$cg),function(site_i){
  
  lapply(unique(df$method_input_code), function(method_input_code_i){
    
    # Subset the data - keeping only one set of method / SNP set
    sub_data <- df %>% filter(method_input_code == method_input_code_i & cg == site_i)
  
    mod <- readRDS(file=here(paste0("outputs/ValidationCommonGarden/HeightModels/stan_models/M3/",
                                    site_i,"_",method_input_code_i,".rds")))[["mod"]]
                     
    # Save coefficients
    broom.mixed::tidyMCMC(mod,
                          droppars = NULL, 
                          robust = FALSE, # give mean and standard deviation
                          ess = F, 
                          rhat = F, 
                          conf.int = T,
                          conf.level = 0.95) %>% 
      #filter(str_detect(term, c('beta'))) %>% 
      dplyr::rename(mean=estimate,
                    std_deviation=std.error,
                    conf_low=conf.low,
                    conf_high=conf.high) %>% 
      mutate(cg = site_i,
             method_input_code = method_input_code_i,
             method_input_name = unique(sub_data$method_input_name),
             method = unique(sub_data$method),
             input_name = unique(sub_data$input_name),
             input_code = unique(sub_data$input_code),
             .before=1)     

}) %>% bind_rows()
 
}) %>% bind_rows()
  
coefftab %>% 
  saveRDS(file=here("outputs/ValidationCommonGarden/HeightModels/coefftab_M3.rds"))

4.3.3 Model coefficients

We plot the mean and 95% credible intervals of the \(\beta_{X_1}\), \(\beta_{X_2}\) and \(\beta_H\) coefficients. Graph titles include the time in months corresponding to the age at which height and survival were recorded. Coefficients in the green area have the expected sign, reflecting lower height in populations with higher GO predictions.

Code
generate_interval_plots(model_i = "M3")
Joining with `by = join_by(cg, method_input_code, method_input_name, method, input_name, input_code)`
Warning in geom_rect(inherit.aes = FALSE, aes(xmin = -Inf, xmax = Inf, ymin = -Inf, : All aesthetics have length 1, but the data has 185 rows.
ℹ Please consider using `annotate()` or provide this layer with data containing a single row.
Warning in geom_rect(inherit.aes = FALSE, aes(xmin = -Inf, xmax = Inf, ymin = -Inf, : All aesthetics have length 1, but the data has 150 rows.
ℹ Please consider using `annotate()` or provide this layer with data containing a single row.
Joining with `by = join_by(cg, method_input_code, method_input_name, method, input_name, input_code)`
Joining with `by = join_by(cg, method_input_code, method_input_name, method, input_name, input_code)`
Joining with `by = join_by(cg, method_input_code, method_input_name, method, input_name, input_code)`
[[1]]
Warning in geom_rect(inherit.aes = FALSE, aes(xmin = -Inf, xmax = Inf, ymin = -Inf, : All aesthetics have length 1, but the data has 185 rows.
ℹ Please consider using `annotate()` or provide this layer with data containing a single row.


[[2]]


[[3]]


[[4]]

4.3.4 Predicted quadratic relationships

We plot the predicted quadratic relationships between tree height and GO predictions / CTD.

Code
generate_poly_plots(model_i = "M3")

4.4 Model 4

4.4.1 Model equation and code

Model 4 = Model 2 + Initial height as a confounder + Clone varying intercepts

\[\begin{align*} Y_{ipb} &\sim \mathcal{N}(\mu_{pb},\sigma^2) \\ \mu_{pb} &= \beta_0 + B_b + \beta_{X1}X_p + \beta_{X2}X^2_p + \beta_H H_p\\ \beta_0 &\sim \mathcal{N}(\text{mean}(Y),2) \\ C_c &\sim \mathcal{N}(0,\sigma_C^2) \\ \begin{bmatrix} \sigma \\ \sigma_C \end{bmatrix} &\sim \text{Exp}(1) \\ \begin{bmatrix} B_b \\ \beta_{X1} \\ \beta_{X2} \\ \beta_H \end{bmatrix} &\sim \mathcal{N}(0,1) \\ \end{align*}\]

  • \(Y_{ipb}\) is the height of the individual \(i\) in the population \(p\) and block \(b\).
  • \(\sigma^{2}\) is the residual variance of the model.
  • \(B_b\) are the block intercepts.
  • \(X_p\) is the GO or CTD of the population \(p\), with \(\beta_{X1}\) and \(\beta_{X2}\) being its linear and quadratic coefficients, respectively. The quadratic term was included to allow for potential nonlinearity in the response, following Fitzpatrick et al. (2021).
  • \(H_p\) is a proxy of the initial height of the trees from population \(p\), i.e. when trees were planted. For that, we used the population varying intercepts calculated across all common gardens in the model 1 of Archambeau et al. (2022).
  • \(\beta_0\) is the model global intercept.
  • \(C_c\) are the clone varying intercepts which follow a normal distribution of variance \(\sigma_C^2\).
Code
stancode = stan_model(here("scripts/StanModels/ValidationCommonGarden_HeightModel_M4.stan"))
print(stancode)
S4 class stanmodel 'anon_model' coded as follows:
data {
  int N;
  vector[N] Y;                              // Response variable (individual height)
  vector[N] X;                              // Genomic offset or climatic transfer distances
  vector[N] H;                              // Initial height of the populations
  int<lower=0> nb_bloc;                     // Number of blocks
  int<lower=0, upper=nb_bloc> bloc[N];      // Blocks
  int<lower=0> nb_clon;                     // Number of clones
  int<lower=0, upper=nb_clon> clon[N];      // Clones
}
parameters {
  vector[nb_bloc] alpha_bloc;               // Intercepts of the blocks
  real beta_0;                              // Global intercept
  real beta_X1;                             // Linear coefficent of GO or CTD
  real beta_X2;                             // Quadratic coefficent of GO or CTD
  real beta_H;                              // Coefficient of the initial height of the populations
  real<lower = 0>  sigma;                   // Residual variance of the model
  real<lower = 0> sigma_clon;               // Variance of the clone intercepts
  vector[nb_clon] z_clon;                     // Vector for non-centered parameterization

}
transformed parameters {
  vector[N] mu;                             // Linear predictor
  real R_squared;                           // R^2 to evaluate the goodness of fit of the model
  vector[nb_clon] alpha_clon;               // Intercepts of the clones

  alpha_clon = z_clon*sigma_clon;
  mu = beta_0 + alpha_bloc[bloc] + alpha_clon[clon] + beta_X1 * X + beta_X2 * square(X) + beta_H * H;
  R_squared = 1 - variance(Y - mu) / variance(Y);
}
model {

  // Likelihood
  Y ~ normal(mu, sigma);

  // Priors
  beta_0 ~ normal(mean(Y),2);
  sigma ~ exponential(1);
  sigma_clon ~ exponential(1);
  alpha_bloc ~ std_normal();
  z_clon ~ std_normal();
  beta_H ~ std_normal();
  beta_X1 ~ std_normal();
  beta_X2 ~ std_normal();
}
// generated quantities{
//   // log likelihood for loo
//   vector[N] log_lik;
//   vector[N] muhat;
//   for (n in 1:N) {
//     log_lik[n] = normal_lpdf(Y[n] |mu[n],sigma);
//     muhat[n] = normal_rng(mu[n], sigma);
//   }
// } 
Code
params_to_estimate <- c("beta_X1","beta_X2","beta_H","R_squared","sigma")
Code
lapply(unique(height_data$cg), function(site_i){
  
  lapply(unique(df$method_input_code), function(method_input_code_i){
  
    # Subset the data
    df_sub <- df %>% filter(method_input_code == method_input_code_i & cg == site_i)
    
    sub_data <- height_data %>% 
      filter(cg == site_i) %>% 
      left_join(df_sub, by = c("cg","pop")) %>% 
      left_join(pop_heights %>% dplyr::select(any_of(c("height", "pop"))) %>% dplyr::rename(init_height=height), by="pop")
      
    # Data in a list for Stan
    stanlist <- list(N = nrow(sub_data),
                     Y=(sub_data$height-mean(sub_data$height))/sd(sub_data$height),
                     X=(sub_data$varX -mean(sub_data$varX))/sd(sub_data$varX),
                     H=(sub_data$init_height -mean(sub_data$init_height))/sd(sub_data$init_height),
                     nb_bloc = length(unique(sub_data$block)),
                     nb_clon = length(unique(sub_data$clon)),
                     bloc = as.numeric(as.factor(sub_data$block)),
                     clon = as.numeric(as.factor(sub_data$clon)))

    # Run the models
    mod <- sampling(stancode, 
                    data = stanlist, 
                    iter = 2000, 
                    warmup = 1400, 
                    chains = 4, 
                    cores = 4, 
                    save_warmup = FALSE,
                    pars = params_to_estimate) 
    
    # Save the model and the stanlist
    list(mod = mod, 
         stanlist = stanlist) %>% 
      saveRDS(file=here(paste0("outputs/ValidationCommonGarden/HeightModels/stan_models/M4/",site_i,"_",method_input_code_i,".rds")))
    
    
  })
})

Warning message: “Warning: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable. Running the chains for more iterations may help. See https://mc-stan.org/misc/warnings.html#bulk-ess”

Code
coefftab <- lapply(unique(height_data$cg),function(site_i){
  
  lapply(unique(df$method_input_code), function(method_input_code_i){
    
    # Subset the data - keeping only one set of method / SNP set
    sub_data <- df %>% filter(method_input_code == method_input_code_i & cg == site_i)
  
    mod <- readRDS(file=here(paste0("outputs/ValidationCommonGarden/HeightModels/stan_models/M4/",
                                    site_i,"_",method_input_code_i,".rds")))[["mod"]]
                     
    # Save coefficients
    broom.mixed::tidyMCMC(mod,
                          droppars = NULL, 
                          robust = FALSE, # give mean and standard deviation
                          ess = F, 
                          rhat = F, 
                          conf.int = T,
                          conf.level = 0.95) %>% 
      dplyr::rename(mean=estimate,
                    std_deviation=std.error,
                    conf_low=conf.low,
                    conf_high=conf.high) %>% 
      mutate(cg = site_i,
             method_input_code = method_input_code_i,
             method_input_name = unique(sub_data$method_input_name),
             method = unique(sub_data$method),
             input_name = unique(sub_data$input_name),
             input_code = unique(sub_data$input_code),
             .before=1)     

}) %>% bind_rows()
 
}) %>% bind_rows()
  
coefftab %>% 
  saveRDS(file=here("outputs/ValidationCommonGarden/HeightModels/coefftab_M4.rds"))

4.4.2 Model coefficients

We plot the mean and 95% credible intervals of the \(\beta_{X_1}\), \(\beta_{X_2}\) and \(\beta_H\) coefficients. Graph titles include the time in months corresponding to the age at which height and survival were recorded. Coefficients in the green area have the expected sign, reflecting lower height in populations with higher GO predictions.

Code
generate_interval_plots(model_i = "M4")
Joining with `by = join_by(cg, method_input_code, method_input_name, method, input_name, input_code)`
Warning in geom_rect(inherit.aes = FALSE, aes(xmin = -Inf, xmax = Inf, ymin = -Inf, : All aesthetics have length 1, but the data has 185 rows.
ℹ Please consider using `annotate()` or provide this layer with data containing a single row.
Warning in geom_rect(inherit.aes = FALSE, aes(xmin = -Inf, xmax = Inf, ymin = -Inf, : All aesthetics have length 1, but the data has 150 rows.
ℹ Please consider using `annotate()` or provide this layer with data containing a single row.
Joining with `by = join_by(cg, method_input_code, method_input_name, method, input_name, input_code)`
Joining with `by = join_by(cg, method_input_code, method_input_name, method, input_name, input_code)`
Joining with `by = join_by(cg, method_input_code, method_input_name, method, input_name, input_code)`
[[1]]
Warning in geom_rect(inherit.aes = FALSE, aes(xmin = -Inf, xmax = Inf, ymin = -Inf, : All aesthetics have length 1, but the data has 185 rows.
ℹ Please consider using `annotate()` or provide this layer with data containing a single row.


[[2]]


[[3]]


[[4]]

4.4.3 Predicted quadratic relationships

We plot the predicted quadratic relationships between tree height and GO predictions / CTD.

Code
generate_poly_plots(model_i = "M4")

4.5 Sample size for each clone

We look at the number of height measurements for each clone. In the tables below (one for each common garden), the first column correspond to the number of height measurements and the second column correspond to the number of clones with this number of measurements.

Code
lapply(unique(height_data$cg), function(site_i){
  
  
height_data %>% 
    filter(cg==site_i) %>% 
    dplyr::select(pop,clon) %>% 
    drop_na() %>% 
    group_by(clon) %>% 
    dplyr::summarise(nb_measurements=n()) %>% 
    count(nb_measurements)
  
  
}) %>% set_names(unique(height_data$cg))
$asturias
# A tibble: 5 × 2
  nb_measurements     n
            <int> <int>
1               4     1
2               5     5
3               6    30
4               7   129
5               8   358

$bordeaux
# A tibble: 5 × 2
  nb_measurements     n
            <int> <int>
1               4     3
2               5    13
3               6    36
4               7   108
5               8   271

$caceres
# A tibble: 4 × 2
  nb_measurements     n
            <int> <int>
1               1   183
2               2    58
3               3    11
4               4     2

$madrid
# A tibble: 6 × 2
  nb_measurements     n
            <int> <int>
1               1   148
2               2   149
3               3   101
4               4    51
5               5    15
6               6     3

$portugal
# A tibble: 8 × 2
  nb_measurements     n
            <int> <int>
1               1     7
2               2    20
3               3    49
4               4    91
5               5   135
6               6   127
7               7    66
8               8    26

5 DRYAD repository

DRYAD REPOSITORY: We export mortality and height data in the DRYAD repository: CommonGardenData_cleaned.csv dataset.

Code
pheno_data %>% 
  dplyr::rename(cg = site) %>% 
  dplyr::select(cg,block,pop,clon,tree,any_of(height_measurements),any_of(surv_measurements)) %>% 
  write_csv(here("data/DryadRepo/CommonGardenData_cleaned.csv"))

6 Session information

Code
devtools::session_info()
─ Session info ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
 setting  value
 version  R version 4.4.2 (2024-10-31)
 os       Ubuntu 22.04.5 LTS
 system   x86_64, linux-gnu
 ui       X11
 language (EN)
 collate  fr_FR.UTF-8
 ctype    fr_FR.UTF-8
 tz       Europe/Paris
 date     2025-01-16
 pandoc   3.1.11 @ /usr/lib/rstudio/resources/app/bin/quarto/bin/tools/x86_64/ (via rmarkdown)

─ Packages ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
 package        * version   date (UTC) lib source
 abind            1.4-5     2016-07-21 [1] CRAN (R 4.4.0)
 arrayhelpers     1.1-0     2020-02-04 [1] CRAN (R 4.4.0)
 backports        1.5.0     2024-05-23 [1] CRAN (R 4.4.0)
 bayesplot      * 1.11.1    2024-02-15 [1] CRAN (R 4.4.0)
 bit              4.0.5     2022-11-15 [1] CRAN (R 4.4.0)
 bit64            4.0.5     2020-08-30 [1] CRAN (R 4.4.0)
 bitops           1.0-7     2021-04-24 [1] CRAN (R 4.4.0)
 broom            1.0.6     2024-05-17 [1] CRAN (R 4.4.0)
 broom.mixed      0.2.9.5   2024-04-01 [1] CRAN (R 4.4.0)
 cachem           1.0.8     2023-05-01 [1] CRAN (R 4.4.0)
 car              3.1-2     2023-03-30 [1] CRAN (R 4.4.0)
 carData          3.0-5     2022-01-06 [1] CRAN (R 4.4.0)
 cellranger       1.1.0     2016-07-27 [1] CRAN (R 4.4.0)
 checkmate        2.3.1     2023-12-04 [1] CRAN (R 4.4.0)
 class            7.3-23    2025-01-01 [4] CRAN (R 4.4.2)
 classInt         0.4-10    2023-09-05 [1] CRAN (R 4.4.0)
 cli              3.6.2     2023-12-11 [1] CRAN (R 4.4.0)
 coda             0.19-4.1  2024-01-31 [1] CRAN (R 4.4.0)
 codetools        0.2-19    2023-02-01 [4] CRAN (R 4.2.2)
 colorspace       2.1-0     2023-01-23 [1] CRAN (R 4.4.0)
 corrplot         0.92      2021-11-18 [1] CRAN (R 4.4.0)
 cowplot        * 1.1.3     2024-01-22 [1] CRAN (R 4.4.0)
 crayon           1.5.2     2022-09-29 [1] CRAN (R 4.4.0)
 DBI              1.2.2     2024-02-16 [1] CRAN (R 4.4.0)
 devtools         2.4.5     2022-10-11 [1] CRAN (R 4.4.0)
 digest           0.6.35    2024-03-11 [1] CRAN (R 4.4.0)
 distributional   0.4.0     2024-02-07 [1] CRAN (R 4.4.0)
 dplyr          * 1.1.4     2023-11-17 [1] CRAN (R 4.4.0)
 e1071            1.7-14    2023-12-06 [1] CRAN (R 4.4.0)
 ellipsis         0.3.2     2021-04-29 [1] CRAN (R 4.4.0)
 evaluate         0.23      2023-11-01 [1] CRAN (R 4.4.0)
 fansi            1.0.6     2023-12-08 [1] CRAN (R 4.4.0)
 farver           2.1.2     2024-05-13 [1] CRAN (R 4.4.0)
 fastmap          1.1.1     2023-02-24 [1] CRAN (R 4.4.0)
 forcats        * 1.0.0     2023-01-29 [1] CRAN (R 4.4.0)
 fs               1.6.4     2024-04-25 [1] CRAN (R 4.4.0)
 furrr            0.3.1     2022-08-15 [1] CRAN (R 4.4.0)
 future           1.33.2    2024-03-26 [1] CRAN (R 4.4.0)
 generics         0.1.3     2022-07-05 [1] CRAN (R 4.4.0)
 ggdist           3.3.2     2024-03-05 [1] CRAN (R 4.4.0)
 ggplot2        * 3.5.1     2024-04-23 [1] CRAN (R 4.4.0)
 ggpubr         * 0.6.0.999 2024-06-05 [1] Github (kassambara/ggpubr@6aeb4f7)
 ggsignif         0.6.4     2022-10-13 [1] CRAN (R 4.4.0)
 globals          0.16.3    2024-03-08 [1] CRAN (R 4.4.0)
 glue             1.7.0     2024-01-09 [1] CRAN (R 4.4.0)
 gridExtra        2.3       2017-09-09 [1] CRAN (R 4.4.0)
 gtable           0.3.5     2024-04-22 [1] CRAN (R 4.4.0)
 HDInterval       0.2.4     2022-11-17 [1] CRAN (R 4.4.0)
 here           * 1.0.1     2020-12-13 [1] CRAN (R 4.4.0)
 highr            0.10      2022-12-22 [1] CRAN (R 4.4.0)
 hms              1.1.3     2023-03-21 [1] CRAN (R 4.4.0)
 htmltools        0.5.8.1   2024-04-04 [1] CRAN (R 4.4.0)
 htmlwidgets      1.6.4     2023-12-06 [1] CRAN (R 4.4.0)
 httpuv           1.6.15    2024-03-26 [1] CRAN (R 4.4.0)
 httr             1.4.7     2023-08-15 [1] CRAN (R 4.4.0)
 inline           0.3.19    2021-05-31 [1] CRAN (R 4.4.0)
 jsonlite         1.8.8     2023-12-04 [1] CRAN (R 4.4.0)
 kableExtra     * 1.4.0     2024-01-24 [1] CRAN (R 4.4.0)
 KernSmooth       2.23-26   2025-01-01 [4] CRAN (R 4.4.2)
 knitr          * 1.46      2024-04-06 [1] CRAN (R 4.4.0)
 labeling         0.4.3     2023-08-29 [1] CRAN (R 4.4.0)
 later            1.3.2     2023-12-06 [1] CRAN (R 4.4.0)
 latex2exp      * 0.9.6     2022-11-28 [1] CRAN (R 4.4.0)
 lattice          0.22-5    2023-10-24 [4] CRAN (R 4.3.1)
 lifecycle        1.0.4     2023-11-07 [1] CRAN (R 4.4.0)
 listenv          0.9.1     2024-01-29 [1] CRAN (R 4.4.0)
 loo              2.7.0     2024-02-24 [1] CRAN (R 4.4.0)
 lubridate      * 1.9.3     2023-09-27 [1] CRAN (R 4.4.0)
 magrittr       * 2.0.3     2022-03-30 [1] CRAN (R 4.4.0)
 matrixStats      1.3.0     2024-04-11 [1] CRAN (R 4.4.0)
 memoise          2.0.1     2021-11-26 [1] CRAN (R 4.4.0)
 mime             0.12      2021-09-28 [1] CRAN (R 4.4.0)
 miniUI           0.1.1.1   2018-05-18 [1] CRAN (R 4.4.0)
 munsell          0.5.1     2024-04-01 [1] CRAN (R 4.4.0)
 nlme             3.1-166   2024-08-14 [4] CRAN (R 4.4.2)
 paletteer      * 1.6.0     2024-01-21 [1] CRAN (R 4.4.0)
 parallelly       1.37.1    2024-02-29 [1] CRAN (R 4.4.0)
 pillar           1.9.0     2023-03-22 [1] CRAN (R 4.4.0)
 pkgbuild         1.4.4     2024-03-17 [1] CRAN (R 4.4.0)
 pkgconfig        2.0.3     2019-09-22 [1] CRAN (R 4.4.0)
 pkgload          1.3.4     2024-01-16 [1] CRAN (R 4.4.0)
 plyr             1.8.9     2023-10-02 [1] CRAN (R 4.4.0)
 png            * 0.1-8     2022-11-29 [1] CRAN (R 4.4.0)
 posterior        1.5.0     2023-10-31 [1] CRAN (R 4.4.0)
 profvis          0.3.8     2023-05-02 [1] CRAN (R 4.4.0)
 promises         1.3.0     2024-04-05 [1] CRAN (R 4.4.0)
 proxy            0.4-27    2022-06-09 [1] CRAN (R 4.4.0)
 purrr          * 1.0.2     2023-08-10 [1] CRAN (R 4.4.0)
 QuickJSR         1.1.3     2024-01-31 [1] CRAN (R 4.4.0)
 R6               2.5.1     2021-08-19 [1] CRAN (R 4.4.0)
 ragg             1.3.0     2024-03-13 [1] CRAN (R 4.4.0)
 raster         * 3.6-26    2023-10-14 [1] CRAN (R 4.4.0)
 RColorBrewer   * 1.1-3     2022-04-03 [1] CRAN (R 4.4.0)
 Rcpp             1.0.12    2024-01-09 [1] CRAN (R 4.4.0)
 RcppParallel     5.1.7     2023-02-27 [1] CRAN (R 4.4.0)
 RCurl          * 1.98-1.14 2024-01-09 [1] CRAN (R 4.4.0)
 readr          * 2.1.5     2024-01-10 [1] CRAN (R 4.4.0)
 readxl         * 1.4.3     2023-07-06 [1] CRAN (R 4.4.0)
 rematch2         2.1.2     2020-05-01 [1] CRAN (R 4.4.0)
 remotes          2.5.0     2024-03-17 [1] CRAN (R 4.4.0)
 reshape2       * 1.4.4     2020-04-09 [1] CRAN (R 4.4.0)
 rlang            1.1.4     2024-06-04 [1] CRAN (R 4.4.0)
 rmarkdown        2.26      2024-03-05 [1] CRAN (R 4.4.0)
 rnaturalearth  * 1.0.1     2023-12-15 [1] CRAN (R 4.4.0)
 rprojroot        2.0.4     2023-11-05 [1] CRAN (R 4.4.0)
 rstan          * 2.32.6    2024-03-05 [1] CRAN (R 4.4.0)
 rstatix          0.7.2     2023-02-01 [1] CRAN (R 4.4.0)
 rstudioapi       0.16.0    2024-03-24 [1] CRAN (R 4.4.0)
 scales           1.3.0     2023-11-28 [1] CRAN (R 4.4.0)
 sessioninfo      1.2.2     2021-12-06 [1] CRAN (R 4.4.0)
 sf               1.0-16    2024-03-24 [1] CRAN (R 4.4.0)
 shiny            1.8.1.1   2024-04-02 [1] CRAN (R 4.4.0)
 sp             * 2.1-4     2024-04-30 [1] CRAN (R 4.4.0)
 StanHeaders    * 2.32.7    2024-04-25 [1] CRAN (R 4.4.0)
 stringi          1.8.4     2024-05-06 [1] CRAN (R 4.4.0)
 stringr        * 1.5.1     2023-11-14 [1] CRAN (R 4.4.0)
 svglite          2.1.3     2023-12-08 [1] CRAN (R 4.4.0)
 svUnit           1.0.6     2021-04-19 [1] CRAN (R 4.4.0)
 systemfonts      1.0.6     2024-03-07 [1] CRAN (R 4.4.0)
 tensorA          0.36.2.1  2023-12-13 [1] CRAN (R 4.4.0)
 terra            1.7-71    2024-01-31 [1] CRAN (R 4.4.0)
 textshaping      0.3.7     2023-10-09 [1] CRAN (R 4.4.0)
 tibble         * 3.2.1     2023-03-20 [1] CRAN (R 4.4.0)
 tidybayes      * 3.0.6     2023-08-12 [1] CRAN (R 4.4.0)
 tidyr          * 1.3.1     2024-01-24 [1] CRAN (R 4.4.0)
 tidyselect       1.2.1     2024-03-11 [1] CRAN (R 4.4.0)
 tidyverse      * 2.0.0     2023-02-22 [1] CRAN (R 4.4.0)
 timechange       0.3.0     2024-01-18 [1] CRAN (R 4.4.0)
 tzdb             0.4.0     2023-05-12 [1] CRAN (R 4.4.0)
 units            0.8-5     2023-11-28 [1] CRAN (R 4.4.0)
 urlchecker       1.0.1     2021-11-30 [1] CRAN (R 4.4.0)
 usethis          2.2.3     2024-02-19 [1] CRAN (R 4.4.0)
 utf8             1.2.4     2023-10-22 [1] CRAN (R 4.4.0)
 vctrs            0.6.5     2023-12-01 [1] CRAN (R 4.4.0)
 viridisLite      0.4.2     2023-05-02 [1] CRAN (R 4.4.0)
 vroom            1.6.5     2023-12-05 [1] CRAN (R 4.4.0)
 withr            3.0.0     2024-01-16 [1] CRAN (R 4.4.0)
 xfun             0.43      2024-03-25 [1] CRAN (R 4.4.0)
 xml2             1.3.6     2023-12-04 [1] CRAN (R 4.4.0)
 xtable         * 1.8-4     2019-04-21 [1] CRAN (R 4.4.0)
 yaml             2.3.8     2023-12-11 [1] CRAN (R 4.4.0)

 [1] /home/juliette/R/x86_64-pc-linux-gnu-library/4.4
 [2] /usr/local/lib/R/site-library
 [3] /usr/lib/R/site-library
 [4] /usr/lib/R/library

──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────

References

Archambeau, Juliette, Marta Benito Garzón, Frédéric Barraquand, Marina de Miguel, Christophe Plomion, and Santiago C González-Martı́nez. 2022. “Combining Climatic and Genomic Data Improves Range-Wide Tree Height Growth Prediction in a Forest Tree.” The American Naturalist 200 (4): E141–59. https://www.journals.uchicago.edu/doi/abs/10.1086/720619.
Fitzpatrick, Matthew C, Vikram E Chhatre, Raju Y Soolanayakanahally, and Stephen R Keller. 2021. “Experimental Support for Genomic Prediction of Climate Maladaptation Using the Machine Learning Approach Gradient Forests.” Molecular Ecology Resources. https://onlinelibrary.wiley.com/doi/abs/10.1111/1755-0998.13374.